Updating for version 0.7.0

Datasets has changed from a single number to arrays

.github/workflows/deployment.yaml

@@ -0,0 +1,25 @@
+name: Deployment
+
+on:
+  push:
+    tags:
+      - '*'
+
+jobs:
+  publish-conda-package:
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v2
+
+    - name: Prepare
+      run: |
+        $CONDA/bin/conda install --quiet --yes conda-build anaconda-client
+        $CONDA/bin/conda config --append channels conda-forge
+        $CONDA/bin/conda config --set anaconda_upload yes
+
+    - name: Build and upload
+      env:
+        ANACONDA_TOKEN: ${{ secrets.ANACONDA_TOKEN }}
+      run: |
+        $CONDA/bin/conda build --token $ANACONDA_TOKEN conda-recipe

.travis.yml

@@ -1,31 +0,0 @@
-language: python
-python:
-  - 3.6
-
-# Build only tagged commits
-if: tag IS present
-
-before_install:
-  - wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh
-  - bash miniconda.sh -b -p $HOME/miniconda
-  - export PATH="$HOME/miniconda/bin:$PATH"
-  - conda config --append channels conda-forge
-  - conda config --set always_yes yes
-  - conda config --set anaconda_upload no
-
-install:
-  - conda update -q conda
-  - conda install -q python=$TRAVIS_PYTHON_VERSION conda-build anaconda-client
-
-script:
-  - conda build conda-recipe
-
-deploy:
-  provider: script
-  script: anaconda -t $ANACONDA_TOKEN upload $HOME/miniconda/conda-bld/**/pyzebra-*.tar.bz2
-  on:
-    branch: master
-    tags: true
-
-notifications:
-  email: false

.vscode/launch.json

@@ -8,6 +8,7 @@
             "program": "${workspaceFolder}/pyzebra/app/cli.py",
             "console": "internalConsole",
             "env": {},
+            "justMyCode": false,
         },
     ]
 }

conda-recipe/meta.yaml

@@ -15,19 +15,16 @@ build:
 
 requirements:
   build:
-    - python >=3.6
+    - python >=3.7
     - setuptools
   run:
-    - python >=3.6
+    - python >=3.7
     - numpy
     - scipy
-    - pandas
     - h5py
-    - bokeh
-    - matplotlib
+    - bokeh =2.4
     - numba
-    - lmfit
-    - uncertainties
+    - lmfit >=1.0.2
 
 
 about:

make_release.py

@@ -3,9 +3,15 @@
 import argparse
 import os
 import re
+import subprocess
 
 
 def main():
+    branch = subprocess.check_output("git rev-parse --abbrev-ref HEAD", shell=True).decode().strip()
+    if branch != "master":
+        print("Aborting, not on 'master' branch.")
+        return
+
     filepath = "pyzebra/__init__.py"
 
     parser = argparse.ArgumentParser()

pyzebra/__init__.py

@@ -1,9 +1,9 @@
 from pyzebra.anatric import *
-from pyzebra.ccl_findpeaks import ccl_findpeaks
-from pyzebra.fit2 import fitccl
+from pyzebra.ccl_io import *
+from pyzebra.ccl_process import *
 from pyzebra.h5 import *
-from pyzebra.ccl_io import load_1D, parse_1D, export_comm
-from pyzebra.param_study_moduls import add_dict, auto, merge, scan_dict
+from pyzebra.utils import *
 from pyzebra.xtal import *
+from pyzebra.sxtal_refgen import *
 
-__version__ = "0.1.3"
+__version__ = "0.7.0"

pyzebra/anatric.py

@@ -7,6 +7,7 @@ DATA_FACTORY_IMPLEMENTATION = [
     "morph",
     "d10",
 ]
+
 REFLECTION_PRINTER_FORMATS = [
     "rafin",
     "rafinf",
@@ -20,11 +21,21 @@ REFLECTION_PRINTER_FORMATS = [
     "oksana",
 ]
 
+ANATRIC_PATH = "/afs/psi.ch/project/sinq/rhel7/bin/anatric"
 ALGORITHMS = ["adaptivemaxcog", "adaptivedynamic"]
 
 
-def anatric(config_file, anatric_path="/afs/psi.ch/project/sinq/rhel7/bin/anatric"):
-    subprocess.run([anatric_path, config_file], check=True)
+def anatric(config_file, anatric_path=ANATRIC_PATH, cwd=None):
+    comp_proc = subprocess.run(
+        [anatric_path, config_file],
+        stdout=subprocess.PIPE,
+        stderr=subprocess.STDOUT,
+        cwd=cwd,
+        check=True,
+        text=True,
+    )
+    print(" ".join(comp_proc.args))
+    print(comp_proc.stdout)
 
 
 class AnatricConfig:
@@ -51,10 +62,13 @@ class AnatricConfig:
     def save_as(self, filename):
         self._tree.write(filename)
 
+    def tostring(self):
+        return ET.tostring(self._tree.getroot(), encoding="unicode")
+
     def _get_attr(self, name, tag, attr):
         elem = self._tree.find(name).find(tag)
         if elem is None:
-            return None
+            return ""
         return elem.attrib[attr]
 
     def _set_attr(self, name, tag, attr, value):
@@ -217,7 +231,7 @@ class AnatricConfig:
         elem = self._tree.find("crystal").find("UB")
         if elem is not None:
             return elem.text
-        return None
+        return ""
 
     @crystal_UB.setter
     def crystal_UB(self, value):
@@ -236,12 +250,37 @@ class AnatricConfig:
 
     @property
     def dataFactory_dist1(self):
-        return self._tree.find("DataFactory").find("dist1").attrib["value"]
+        elem = self._tree.find("DataFactory").find("dist1")
+        if elem is not None:
+            return elem.attrib["value"]
+        return ""
 
     @dataFactory_dist1.setter
     def dataFactory_dist1(self, value):
         self._tree.find("DataFactory").find("dist1").attrib["value"] = value
 
+    @property
+    def dataFactory_dist2(self):
+        elem = self._tree.find("DataFactory").find("dist2")
+        if elem is not None:
+            return elem.attrib["value"]
+        return ""
+
+    @dataFactory_dist2.setter
+    def dataFactory_dist2(self, value):
+        self._tree.find("DataFactory").find("dist2").attrib["value"] = value
+
+    @property
+    def dataFactory_dist3(self):
+        elem = self._tree.find("DataFactory").find("dist3")
+        if elem is not None:
+            return elem.attrib["value"]
+        return ""
+
+    @dataFactory_dist3.setter
+    def dataFactory_dist3(self, value):
+        self._tree.find("DataFactory").find("dist3").attrib["value"] = value
+
     @property
     def reflectionPrinter_format(self):
         return self._tree.find("ReflectionPrinter").attrib["format"]
@@ -253,6 +292,14 @@ class AnatricConfig:
 
         self._tree.find("ReflectionPrinter").attrib["format"] = value
 
+    @property
+    def reflectionPrinter_file(self):
+        return self._tree.find("ReflectionPrinter").attrib["file"]
+
+    @reflectionPrinter_file.setter
+    def reflectionPrinter_file(self, value):
+        self._tree.find("ReflectionPrinter").attrib["file"] = value
+
     @property
     def algorithm(self):
         return self._tree.find("Algorithm").attrib["implementation"]
@@ -269,7 +316,7 @@ class AnatricConfig:
     def _get_alg_attr(self, alg, tag, attr):
         param_elem = self._alg_elems[alg].find(tag)
         if param_elem is None:
-            return None
+            return ""
         return param_elem.attrib[attr]
 
     def _set_alg_attr(self, alg, tag, attr, value):

pyzebra/app/app.py

@@ -2,14 +2,19 @@ import logging
 import sys
 from io import StringIO
 
+import pyzebra
 from bokeh.io import curdoc
 from bokeh.layouts import column, row
-from bokeh.models import Tabs, TextAreaInput
+from bokeh.models import Button, Panel, Tabs, TextAreaInput, TextInput
 
 import panel_ccl_integrate
+import panel_ccl_compare
 import panel_hdf_anatric
+import panel_hdf_param_study
 import panel_hdf_viewer
-
+import panel_param_study
+import panel_spind
+import panel_ccl_prepare
 
 doc = curdoc()
 
@@ -23,14 +28,46 @@ bokeh_logger = logging.getLogger("bokeh")
 bokeh_logger.addHandler(bokeh_handler)
 bokeh_log_textareainput = TextAreaInput(title="server output:", height=150)
 
-# Final layout
-tab_hdf_viewer = panel_hdf_viewer.create()
-tab_hdf_anatric = panel_hdf_anatric.create()
-tab_ccl_integrate = panel_ccl_integrate.create()
+def proposal_textinput_callback(_attr, _old, _new):
+    apply_button.disabled = False
 
+proposal_textinput = TextInput(title="Proposal number:", name="")
+proposal_textinput.on_change("value_input", proposal_textinput_callback)
+doc.proposal_textinput = proposal_textinput
+
+def apply_button_callback():
+    proposal = proposal_textinput.value.strip()
+    if proposal:
+        try:
+            proposal_path = pyzebra.find_proposal_path(proposal)
+        except ValueError as e:
+            print(e)
+            return
+        apply_button.disabled = True
+    else:
+        proposal_path = ""
+
+    proposal_textinput.name = proposal_path
+
+apply_button = Button(label="Apply", button_type="primary")
+apply_button.on_click(apply_button_callback)
+
+# Final layout
 doc.add_root(
     column(
-        Tabs(tabs=[tab_hdf_viewer, tab_hdf_anatric, tab_ccl_integrate]),
+        Tabs(
+            tabs=[
+                Panel(child=column(proposal_textinput, apply_button), title="user config"),
+                panel_hdf_viewer.create(),
+                panel_hdf_anatric.create(),
+                panel_ccl_prepare.create(),
+                panel_ccl_integrate.create(),
+                panel_ccl_compare.create(),
+                panel_param_study.create(),
+                panel_hdf_param_study.create(),
+                panel_spind.create(),
+            ]
+        ),
         row(stdout_textareainput, bokeh_log_textareainput, sizing_mode="scale_both"),
     )
 )

pyzebra/app/cli.py

@@ -6,6 +6,7 @@ from bokeh.application.application import Application
 from bokeh.application.handlers import ScriptHandler
 from bokeh.server.server import Server
 
+from pyzebra import ANATRIC_PATH, SXTAL_REFGEN_PATH
 from pyzebra.app.handler import PyzebraHandler
 
 logging.basicConfig(format="%(asctime)s %(message)s", level=logging.INFO)
@@ -38,10 +39,18 @@ def main():
     )
 
     parser.add_argument(
-        "--anatric-path",
+        "--anatric-path", type=str, default=ANATRIC_PATH, help="path to anatric executable",
+    )
+
+    parser.add_argument(
+        "--sxtal-refgen-path",
         type=str,
-        default="/afs/psi.ch/project/sinq/rhel7/bin/anatric",
-        help="path to anatric executable",
+        default=SXTAL_REFGEN_PATH,
+        help="path to Sxtal_Refgen executable",
+    )
+
+    parser.add_argument(
+        "--spind-path", type=str, default=None, help="path to spind scripts folder",
     )
 
     parser.add_argument(
@@ -55,7 +64,7 @@ def main():
 
     logger.info(app_path)
 
-    pyzebra_handler = PyzebraHandler(args.anatric_path)
+    pyzebra_handler = PyzebraHandler(args.anatric_path, args.spind_path)
     handler = ScriptHandler(filename=app_path, argv=args.args)
     server = Server(
         {"/": Application(pyzebra_handler, handler)},

pyzebra/app/handler.py

@@ -5,7 +5,7 @@ class PyzebraHandler(Handler):
     """Provides a mechanism for generic bokeh applications to build up new streamvis documents.
     """
 
-    def __init__(self, anatric_path):
+    def __init__(self, anatric_path, spind_path):
         """Initialize a pyzebra handler for bokeh applications.
 
         Args:
@@ -14,6 +14,7 @@ class PyzebraHandler(Handler):
         super().__init__()  # no-op
 
         self.anatric_path = anatric_path
+        self.spind_path = spind_path
 
     def modify_document(self, doc):
         """Modify an application document with pyzebra specific features.
@@ -26,5 +27,6 @@ class PyzebraHandler(Handler):
         """
         doc.title = "pyzebra"
         doc.anatric_path = self.anatric_path
+        doc.spind_path = self.spind_path
 
         return doc

pyzebra/app/panel_ccl_compare.py

@@ -0,0 +1,718 @@
+import base64
+import io
+import os
+import tempfile
+import types
+
+import numpy as np
+from bokeh.io import curdoc
+from bokeh.layouts import column, row
+from bokeh.models import (
+    BasicTicker,
+    Button,
+    CellEditor,
+    CheckboxEditor,
+    CheckboxGroup,
+    ColumnDataSource,
+    CustomJS,
+    DataRange1d,
+    DataTable,
+    Div,
+    Dropdown,
+    FileInput,
+    Grid,
+    Legend,
+    Line,
+    LinearAxis,
+    MultiLine,
+    MultiSelect,
+    NumberEditor,
+    Panel,
+    PanTool,
+    Plot,
+    RadioGroup,
+    ResetTool,
+    Scatter,
+    Select,
+    Spacer,
+    Span,
+    Spinner,
+    TableColumn,
+    TextAreaInput,
+    WheelZoomTool,
+    Whisker,
+)
+
+import pyzebra
+from pyzebra.ccl_io import EXPORT_TARGETS
+from pyzebra.ccl_process import AREA_METHODS
+
+
+javaScript = """
+let j = 0;
+for (let i = 0; i < js_data.data['fname'].length; i++) {
+    if (js_data.data['content'][i] === "") continue;
+
+    setTimeout(function() {
+        const blob = new Blob([js_data.data['content'][i]], {type: 'text/plain'})
+        const link = document.createElement('a');
+        document.body.appendChild(link);
+        const url = window.URL.createObjectURL(blob);
+        link.href = url;
+        link.download = js_data.data['fname'][i] + js_data.data['ext'][i];
+        link.click();
+        window.URL.revokeObjectURL(url);
+        document.body.removeChild(link);
+    }, 100 * j)
+
+    j++;
+}
+"""
+
+
+def create():
+    doc = curdoc()
+    dataset1 = []
+    dataset2 = []
+    fit_params = {}
+    js_data = ColumnDataSource(data=dict(content=["", ""], fname=["", ""], ext=["", ""]))
+
+    def file_select_update_for_proposal():
+        proposal_path = proposal_textinput.name
+        if proposal_path:
+            file_list = []
+            for file in os.listdir(proposal_path):
+                if file.endswith((".ccl")):
+                    file_list.append((os.path.join(proposal_path, file), file))
+            file_select.options = file_list
+            file_open_button.disabled = False
+        else:
+            file_select.options = []
+            file_open_button.disabled = True
+
+    doc.add_periodic_callback(file_select_update_for_proposal, 5000)
+
+    def proposal_textinput_callback(_attr, _old, _new):
+        file_select_update_for_proposal()
+
+    proposal_textinput = doc.proposal_textinput
+    proposal_textinput.on_change("name", proposal_textinput_callback)
+
+    def _init_datatable():
+        # dataset2 should have the same metadata as dataset1
+        scan_list = [s["idx"] for s in dataset1]
+        hkl = [f'{s["h"]} {s["k"]} {s["l"]}' for s in dataset1]
+        export = [s["export"] for s in dataset1]
+
+        twotheta = [np.median(s["twotheta"]) if "twotheta" in s else None for s in dataset1]
+        gamma = [np.median(s["gamma"]) if "gamma" in s else None for s in dataset1]
+        omega = [np.median(s["omega"]) if "omega" in s else None for s in dataset1]
+        chi = [np.median(s["chi"]) if "chi" in s else None for s in dataset1]
+        phi = [np.median(s["phi"]) if "phi" in s else None for s in dataset1]
+        nu = [np.median(s["nu"]) if "nu" in s else None for s in dataset1]
+
+        scan_table_source.data.update(
+            scan=scan_list,
+            hkl=hkl,
+            fit=[0] * len(scan_list),
+            export=export,
+            twotheta=twotheta,
+            gamma=gamma,
+            omega=omega,
+            chi=chi,
+            phi=phi,
+            nu=nu,
+        )
+        scan_table_source.selected.indices = []
+        scan_table_source.selected.indices = [0]
+
+        merge_options = [(str(i), f"{i} ({idx})") for i, idx in enumerate(scan_list)]
+        merge_from_select.options = merge_options
+        merge_from_select.value = merge_options[0][0]
+
+    file_select = MultiSelect(title="Select 2 .ccl files:", width=210, height=250)
+
+    def file_open_button_callback():
+        if len(file_select.value) != 2:
+            print("WARNING: Select exactly 2 .ccl files.")
+            return
+
+        new_data1 = []
+        new_data2 = []
+        for ind, f_path in enumerate(file_select.value):
+            with open(f_path) as file:
+                f_name = os.path.basename(f_path)
+                base, ext = os.path.splitext(f_name)
+                try:
+                    file_data = pyzebra.parse_1D(file, ext)
+                except:
+                    print(f"Error loading {f_name}")
+                    return
+
+            pyzebra.normalize_dataset(file_data, monitor_spinner.value)
+            pyzebra.merge_duplicates(file_data)
+
+            if ind == 0:
+                js_data.data.update(fname=[base, base])
+                new_data1 = file_data
+            else:  # ind = 1
+                new_data2 = file_data
+
+        # ignore extra scans at the end of the longest of the two files
+        min_len = min(len(new_data1), len(new_data2))
+        new_data1 = new_data1[:min_len]
+        new_data2 = new_data2[:min_len]
+
+        nonlocal dataset1, dataset2
+        dataset1 = new_data1
+        dataset2 = new_data2
+        _init_datatable()
+
+    file_open_button = Button(label="Open New", width=100, disabled=True)
+    file_open_button.on_click(file_open_button_callback)
+
+    def upload_button_callback(_attr, _old, _new):
+        if len(upload_button.filename) != 2:
+            print("WARNING: Upload exactly 2 .ccl files.")
+            return
+
+        new_data1 = []
+        new_data2 = []
+        for ind, (f_str, f_name) in enumerate(zip(upload_button.value, upload_button.filename)):
+            with io.StringIO(base64.b64decode(f_str).decode()) as file:
+                base, ext = os.path.splitext(f_name)
+                try:
+                    file_data = pyzebra.parse_1D(file, ext)
+                except:
+                    print(f"Error loading {f_name}")
+                    return
+
+            pyzebra.normalize_dataset(file_data, monitor_spinner.value)
+            pyzebra.merge_duplicates(file_data)
+
+            if ind == 0:
+                js_data.data.update(fname=[base, base])
+                new_data1 = file_data
+            else:  # ind = 1
+                new_data2 = file_data
+
+        # ignore extra scans at the end of the longest of the two files
+        min_len = min(len(new_data1), len(new_data2))
+        new_data1 = new_data1[:min_len]
+        new_data2 = new_data2[:min_len]
+
+        nonlocal dataset1, dataset2
+        dataset1 = new_data1
+        dataset2 = new_data2
+        _init_datatable()
+
+    upload_div = Div(text="or upload 2 .ccl files:", margin=(5, 5, 0, 5))
+    upload_button = FileInput(accept=".ccl", multiple=True, width=200)
+    # for on_change("value", ...) or on_change("filename", ...),
+    # see https://github.com/bokeh/bokeh/issues/11461
+    upload_button.on_change("filename", upload_button_callback)
+
+    def monitor_spinner_callback(_attr, old, new):
+        if dataset1 and dataset2:
+            pyzebra.normalize_dataset(dataset1, new)
+            pyzebra.normalize_dataset(dataset2, new)
+            _update_plot()
+
+    monitor_spinner = Spinner(title="Monitor:", mode="int", value=100_000, low=1, width=145)
+    monitor_spinner.on_change("value", monitor_spinner_callback)
+
+    def _update_table():
+        fit_ok = [(1 if "fit" in scan else 0) for scan in dataset1]
+        export = [scan["export"] for scan in dataset1]
+        scan_table_source.data.update(fit=fit_ok, export=export)
+
+    def _update_plot():
+        plot_scatter_source = [plot_scatter1_source, plot_scatter2_source]
+        plot_fit_source = [plot_fit1_source, plot_fit2_source]
+        plot_bkg_source = [plot_bkg1_source, plot_bkg2_source]
+        plot_peak_source = [plot_peak1_source, plot_peak2_source]
+        fit_output = ""
+
+        for ind, scan in enumerate(_get_selected_scan()):
+            scatter_source = plot_scatter_source[ind]
+            fit_source = plot_fit_source[ind]
+            bkg_source = plot_bkg_source[ind]
+            peak_source = plot_peak_source[ind]
+            scan_motor = scan["scan_motor"]
+
+            y = scan["counts"]
+            y_err = scan["counts_err"]
+            x = scan[scan_motor]
+
+            plot.axis[0].axis_label = scan_motor
+            scatter_source.data.update(x=x, y=y, y_upper=y + y_err, y_lower=y - y_err)
+
+            fit = scan.get("fit")
+            if fit is not None:
+                x_fit = np.linspace(x[0], x[-1], 100)
+                fit_source.data.update(x=x_fit, y=fit.eval(x=x_fit))
+
+                x_bkg = []
+                y_bkg = []
+                xs_peak = []
+                ys_peak = []
+                comps = fit.eval_components(x=x_fit)
+                for i, model in enumerate(fit_params):
+                    if "linear" in model:
+                        x_bkg = x_fit
+                        y_bkg = comps[f"f{i}_"]
+
+                    elif any(val in model for val in ("gaussian", "voigt", "pvoigt")):
+                        xs_peak.append(x_fit)
+                        ys_peak.append(comps[f"f{i}_"])
+
+                bkg_source.data.update(x=x_bkg, y=y_bkg)
+                peak_source.data.update(xs=xs_peak, ys=ys_peak)
+                if fit_output:
+                    fit_output = fit_output + "\n\n"
+                fit_output = fit_output + fit.fit_report()
+
+            else:
+                fit_source.data.update(x=[], y=[])
+                bkg_source.data.update(x=[], y=[])
+                peak_source.data.update(xs=[], ys=[])
+
+        fit_output_textinput.value = fit_output
+
+    # Main plot
+    plot = Plot(
+        x_range=DataRange1d(),
+        y_range=DataRange1d(only_visible=True),
+        plot_height=470,
+        plot_width=700,
+    )
+
+    plot.add_layout(LinearAxis(axis_label="Counts"), place="left")
+    plot.add_layout(LinearAxis(axis_label="Scan motor"), place="below")
+
+    plot.add_layout(Grid(dimension=0, ticker=BasicTicker()))
+    plot.add_layout(Grid(dimension=1, ticker=BasicTicker()))
+
+    plot_scatter1_source = ColumnDataSource(dict(x=[0], y=[0], y_upper=[0], y_lower=[0]))
+    plot_scatter1 = plot.add_glyph(
+        plot_scatter1_source, Scatter(x="x", y="y", line_color="steelblue", fill_color="steelblue")
+    )
+    plot.add_layout(
+        Whisker(source=plot_scatter1_source, base="x", upper="y_upper", lower="y_lower")
+    )
+
+    plot_scatter2_source = ColumnDataSource(dict(x=[0], y=[0], y_upper=[0], y_lower=[0]))
+    plot_scatter2 = plot.add_glyph(
+        plot_scatter2_source, Scatter(x="x", y="y", line_color="firebrick", fill_color="firebrick")
+    )
+    plot.add_layout(
+        Whisker(source=plot_scatter2_source, base="x", upper="y_upper", lower="y_lower")
+    )
+
+    plot_fit1_source = ColumnDataSource(dict(x=[0], y=[0]))
+    plot_fit1 = plot.add_glyph(plot_fit1_source, Line(x="x", y="y"))
+
+    plot_fit2_source = ColumnDataSource(dict(x=[0], y=[0]))
+    plot_fit2 = plot.add_glyph(plot_fit2_source, Line(x="x", y="y"))
+
+    plot_bkg1_source = ColumnDataSource(dict(x=[0], y=[0]))
+    plot_bkg1 = plot.add_glyph(
+        plot_bkg1_source, Line(x="x", y="y", line_color="steelblue", line_dash="dashed")
+    )
+
+    plot_bkg2_source = ColumnDataSource(dict(x=[0], y=[0]))
+    plot_bkg2 = plot.add_glyph(
+        plot_bkg2_source, Line(x="x", y="y", line_color="firebrick", line_dash="dashed")
+    )
+
+    plot_peak1_source = ColumnDataSource(dict(xs=[[0]], ys=[[0]]))
+    plot_peak1 = plot.add_glyph(
+        plot_peak1_source, MultiLine(xs="xs", ys="ys", line_color="steelblue", line_dash="dashed")
+    )
+
+    plot_peak2_source = ColumnDataSource(dict(xs=[[0]], ys=[[0]]))
+    plot_peak2 = plot.add_glyph(
+        plot_peak2_source, MultiLine(xs="xs", ys="ys", line_color="firebrick", line_dash="dashed")
+    )
+
+    fit_from_span = Span(location=None, dimension="height", line_dash="dashed")
+    plot.add_layout(fit_from_span)
+
+    fit_to_span = Span(location=None, dimension="height", line_dash="dashed")
+    plot.add_layout(fit_to_span)
+
+    plot.add_layout(
+        Legend(
+            items=[
+                ("data 1", [plot_scatter1]),
+                ("data 2", [plot_scatter2]),
+                ("best fit 1", [plot_fit1]),
+                ("best fit 2", [plot_fit2]),
+                ("peak 1", [plot_peak1]),
+                ("peak 2", [plot_peak2]),
+                ("linear 1", [plot_bkg1]),
+                ("linear 2", [plot_bkg2]),
+            ],
+            location="top_left",
+            click_policy="hide",
+        )
+    )
+
+    plot.add_tools(PanTool(), WheelZoomTool(), ResetTool())
+    plot.toolbar.logo = None
+
+    # Scan select
+    def scan_table_select_callback(_attr, old, new):
+        if not new:
+            # skip empty selections
+            return
+
+        # Avoid selection of multiple indicies (via Shift+Click or Ctrl+Click)
+        if len(new) > 1:
+            # drop selection to the previous one
+            scan_table_source.selected.indices = old
+            return
+
+        if len(old) > 1:
+            # skip unnecessary update caused by selection drop
+            return
+
+        _update_plot()
+
+    def scan_table_source_callback(_attr, _old, new):
+        # unfortunately, we don't know if the change comes from data update or user input
+        # also `old` and `new` are the same for non-scalars
+        for scan1, scan2, export in zip(dataset1, dataset2, new["export"]):
+            scan1["export"] = export
+            scan2["export"] = export
+        _update_preview()
+
+    scan_table_source = ColumnDataSource(
+        dict(
+            scan=[],
+            hkl=[],
+            fit=[],
+            export=[],
+            twotheta=[],
+            gamma=[],
+            omega=[],
+            chi=[],
+            phi=[],
+            nu=[],
+        )
+    )
+    scan_table_source.on_change("data", scan_table_source_callback)
+    scan_table_source.selected.on_change("indices", scan_table_select_callback)
+
+    scan_table = DataTable(
+        source=scan_table_source,
+        columns=[
+            TableColumn(field="scan", title="Scan", editor=CellEditor(), width=50),
+            TableColumn(field="hkl", title="hkl", editor=CellEditor(), width=100),
+            TableColumn(field="fit", title="Fit", editor=CellEditor(), width=50),
+            TableColumn(field="export", title="Export", editor=CheckboxEditor(), width=50),
+            TableColumn(field="twotheta", title="2theta", editor=CellEditor(), width=50),
+            TableColumn(field="gamma", title="gamma", editor=CellEditor(), width=50),
+            TableColumn(field="omega", title="omega", editor=CellEditor(), width=50),
+            TableColumn(field="chi", title="chi", editor=CellEditor(), width=50),
+            TableColumn(field="phi", title="phi", editor=CellEditor(), width=50),
+            TableColumn(field="nu", title="nu", editor=CellEditor(), width=50),
+        ],
+        width=310,  # +60 because of the index column, but excluding twotheta onwards
+        height=350,
+        autosize_mode="none",
+        editable=True,
+    )
+
+    def _get_selected_scan():
+        ind = scan_table_source.selected.indices[0]
+        return dataset1[ind], dataset2[ind]
+
+    merge_from_select = Select(title="scan:", width=145)
+
+    def merge_button_callback():
+        scan_into1, scan_into2 = _get_selected_scan()
+        scan_from1 = dataset1[int(merge_from_select.value)]
+        scan_from2 = dataset2[int(merge_from_select.value)]
+
+        if scan_into1 is scan_from1:
+            print("WARNING: Selected scans for merging are identical")
+            return
+
+        pyzebra.merge_scans(scan_into1, scan_from1)
+        pyzebra.merge_scans(scan_into2, scan_from2)
+        _update_table()
+        _update_plot()
+
+    merge_button = Button(label="Merge into current", width=145)
+    merge_button.on_click(merge_button_callback)
+
+    def restore_button_callback():
+        scan1, scan2 = _get_selected_scan()
+        pyzebra.restore_scan(scan1)
+        pyzebra.restore_scan(scan2)
+        _update_table()
+        _update_plot()
+
+    restore_button = Button(label="Restore scan", width=145)
+    restore_button.on_click(restore_button_callback)
+
+    def fit_from_spinner_callback(_attr, _old, new):
+        fit_from_span.location = new
+
+    fit_from_spinner = Spinner(title="Fit from:", width=145)
+    fit_from_spinner.on_change("value", fit_from_spinner_callback)
+
+    def fit_to_spinner_callback(_attr, _old, new):
+        fit_to_span.location = new
+
+    fit_to_spinner = Spinner(title="to:", width=145)
+    fit_to_spinner.on_change("value", fit_to_spinner_callback)
+
+    def fitparams_add_dropdown_callback(click):
+        # bokeh requires (str, str) for MultiSelect options
+        new_tag = f"{click.item}-{fitparams_select.tags[0]}"
+        fitparams_select.options.append((new_tag, click.item))
+        fit_params[new_tag] = fitparams_factory(click.item)
+        fitparams_select.tags[0] += 1
+
+    fitparams_add_dropdown = Dropdown(
+        label="Add fit function",
+        menu=[
+            ("Linear", "linear"),
+            ("Gaussian", "gaussian"),
+            ("Voigt", "voigt"),
+            ("Pseudo Voigt", "pvoigt"),
+            # ("Pseudo Voigt1", "pseudovoigt1"),
+        ],
+        width=145,
+    )
+    fitparams_add_dropdown.on_click(fitparams_add_dropdown_callback)
+
+    def fitparams_select_callback(_attr, old, new):
+        # Avoid selection of multiple indicies (via Shift+Click or Ctrl+Click)
+        if len(new) > 1:
+            # drop selection to the previous one
+            fitparams_select.value = old
+            return
+
+        if len(old) > 1:
+            # skip unnecessary update caused by selection drop
+            return
+
+        if new:
+            fitparams_table_source.data.update(fit_params[new[0]])
+        else:
+            fitparams_table_source.data.update(dict(param=[], value=[], vary=[], min=[], max=[]))
+
+    fitparams_select = MultiSelect(options=[], height=120, width=145)
+    fitparams_select.tags = [0]
+    fitparams_select.on_change("value", fitparams_select_callback)
+
+    def fitparams_remove_button_callback():
+        if fitparams_select.value:
+            sel_tag = fitparams_select.value[0]
+            del fit_params[sel_tag]
+            for elem in fitparams_select.options:
+                if elem[0] == sel_tag:
+                    fitparams_select.options.remove(elem)
+                    break
+
+            fitparams_select.value = []
+
+    fitparams_remove_button = Button(label="Remove fit function", width=145)
+    fitparams_remove_button.on_click(fitparams_remove_button_callback)
+
+    def fitparams_factory(function):
+        if function == "linear":
+            params = ["slope", "intercept"]
+        elif function == "gaussian":
+            params = ["amplitude", "center", "sigma"]
+        elif function == "voigt":
+            params = ["amplitude", "center", "sigma", "gamma"]
+        elif function == "pvoigt":
+            params = ["amplitude", "center", "sigma", "fraction"]
+        elif function == "pseudovoigt1":
+            params = ["amplitude", "center", "g_sigma", "l_sigma", "fraction"]
+        else:
+            raise ValueError("Unknown fit function")
+
+        n = len(params)
+        fitparams = dict(
+            param=params, value=[None] * n, vary=[True] * n, min=[None] * n, max=[None] * n,
+        )
+
+        if function == "linear":
+            fitparams["value"] = [0, 1]
+            fitparams["vary"] = [False, True]
+            fitparams["min"] = [None, 0]
+
+        elif function == "gaussian":
+            fitparams["min"] = [0, None, None]
+
+        return fitparams
+
+    fitparams_table_source = ColumnDataSource(dict(param=[], value=[], vary=[], min=[], max=[]))
+    fitparams_table = DataTable(
+        source=fitparams_table_source,
+        columns=[
+            TableColumn(field="param", title="Parameter", editor=CellEditor()),
+            TableColumn(field="value", title="Value", editor=NumberEditor()),
+            TableColumn(field="vary", title="Vary", editor=CheckboxEditor()),
+            TableColumn(field="min", title="Min", editor=NumberEditor()),
+            TableColumn(field="max", title="Max", editor=NumberEditor()),
+        ],
+        height=200,
+        width=350,
+        index_position=None,
+        editable=True,
+        auto_edit=True,
+    )
+
+    # start with `background` and `gauss` fit functions added
+    fitparams_add_dropdown_callback(types.SimpleNamespace(item="linear"))
+    fitparams_add_dropdown_callback(types.SimpleNamespace(item="gaussian"))
+    fitparams_select.value = ["gaussian-1"]  # add selection to gauss
+
+    fit_output_textinput = TextAreaInput(title="Fit results:", width=750, height=200)
+
+    def proc_all_button_callback():
+        for scan in [*dataset1, *dataset2]:
+            if scan["export"]:
+                pyzebra.fit_scan(
+                    scan, fit_params, fit_from=fit_from_spinner.value, fit_to=fit_to_spinner.value
+                )
+                pyzebra.get_area(
+                    scan,
+                    area_method=AREA_METHODS[area_method_radiobutton.active],
+                    lorentz=lorentz_checkbox.active,
+                )
+
+        _update_plot()
+        _update_table()
+
+    proc_all_button = Button(label="Process All", button_type="primary", width=145)
+    proc_all_button.on_click(proc_all_button_callback)
+
+    def proc_button_callback():
+        for scan in _get_selected_scan():
+            pyzebra.fit_scan(
+                scan, fit_params, fit_from=fit_from_spinner.value, fit_to=fit_to_spinner.value
+            )
+            pyzebra.get_area(
+                scan,
+                area_method=AREA_METHODS[area_method_radiobutton.active],
+                lorentz=lorentz_checkbox.active,
+            )
+
+        _update_plot()
+        _update_table()
+
+    proc_button = Button(label="Process Current", width=145)
+    proc_button.on_click(proc_button_callback)
+
+    area_method_div = Div(text="Intensity:", margin=(5, 5, 0, 5))
+    area_method_radiobutton = RadioGroup(labels=["Function", "Area"], active=0, width=145)
+
+    intensity_diff_div = Div(text="Intensity difference:", margin=(5, 5, 0, 5))
+    intensity_diff_radiobutton = RadioGroup(
+        labels=["file1 - file2", "file2 - file1"], active=0, width=145
+    )
+
+    lorentz_checkbox = CheckboxGroup(labels=["Lorentz Correction"], width=145, margin=(13, 5, 5, 5))
+
+    export_preview_textinput = TextAreaInput(title="Export file(s) preview:", width=500, height=400)
+
+    def _update_preview():
+        with tempfile.TemporaryDirectory() as temp_dir:
+            temp_file = temp_dir + "/temp"
+            export_data1 = []
+            export_data2 = []
+            for scan1, scan2 in zip(dataset1, dataset2):
+                if scan1["export"]:
+                    export_data1.append(scan1)
+                    export_data2.append(scan2)
+
+            if intensity_diff_radiobutton.active:
+                export_data1, export_data2 = export_data2, export_data1
+
+            pyzebra.export_ccl_compare(
+                export_data1,
+                export_data2,
+                temp_file,
+                export_target_select.value,
+                hkl_precision=int(hkl_precision_select.value),
+            )
+
+            exported_content = ""
+            file_content = []
+            for ext in EXPORT_TARGETS[export_target_select.value]:
+                fname = temp_file + ext
+                if os.path.isfile(fname):
+                    with open(fname) as f:
+                        content = f.read()
+                        exported_content += f"{ext} file:\n" + content
+                else:
+                    content = ""
+                file_content.append(content)
+
+            js_data.data.update(content=file_content)
+            export_preview_textinput.value = exported_content
+
+    def export_target_select_callback(_attr, _old, new):
+        js_data.data.update(ext=EXPORT_TARGETS[new])
+        _update_preview()
+
+    export_target_select = Select(
+        title="Export target:", options=list(EXPORT_TARGETS.keys()), value="fullprof", width=80
+    )
+    export_target_select.on_change("value", export_target_select_callback)
+    js_data.data.update(ext=EXPORT_TARGETS[export_target_select.value])
+
+    def hkl_precision_select_callback(_attr, _old, _new):
+        _update_preview()
+
+    hkl_precision_select = Select(
+        title="hkl precision:", options=["2", "3", "4"], value="2", width=80
+    )
+    hkl_precision_select.on_change("value", hkl_precision_select_callback)
+
+    save_button = Button(label="Download File(s)", button_type="success", width=200)
+    save_button.js_on_click(CustomJS(args={"js_data": js_data}, code=javaScript))
+
+    fitpeak_controls = row(
+        column(fitparams_add_dropdown, fitparams_select, fitparams_remove_button),
+        fitparams_table,
+        Spacer(width=20),
+        column(
+            fit_from_spinner,
+            lorentz_checkbox,
+            area_method_div,
+            area_method_radiobutton,
+            intensity_diff_div,
+            intensity_diff_radiobutton,
+        ),
+        column(fit_to_spinner, proc_button, proc_all_button),
+    )
+
+    scan_layout = column(
+        scan_table,
+        row(monitor_spinner, column(Spacer(height=19), restore_button)),
+        row(column(Spacer(height=19), merge_button), merge_from_select),
+    )
+
+    import_layout = column(file_select, file_open_button, upload_div, upload_button)
+
+    export_layout = column(
+        export_preview_textinput,
+        row(
+            export_target_select, hkl_precision_select, column(Spacer(height=19), row(save_button))
+        ),
+    )
+
+    tab_layout = column(
+        row(import_layout, scan_layout, plot, Spacer(width=30), export_layout),
+        row(fitpeak_controls, fit_output_textinput),
+    )
+
+    return Panel(child=tab_layout, title="ccl compare")

pyzebra/app/panel_ccl_prepare.py

@@ -0,0 +1,745 @@
+import base64
+import io
+import os
+import subprocess
+import tempfile
+
+import numpy as np
+from bokeh.layouts import column, row
+from bokeh.models import (
+    Arrow,
+    BoxZoomTool,
+    Button,
+    CheckboxGroup,
+    ColumnDataSource,
+    CustomJS,
+    Div,
+    Ellipse,
+    FileInput,
+    Legend,
+    LegendItem,
+    LinearAxis,
+    MultiLine,
+    MultiSelect,
+    NormalHead,
+    NumericInput,
+    Panel,
+    PanTool,
+    Plot,
+    RadioGroup,
+    Range1d,
+    ResetTool,
+    Scatter,
+    Select,
+    Spacer,
+    Spinner,
+    Text,
+    TextAreaInput,
+    TextInput,
+    WheelZoomTool,
+)
+from bokeh.palettes import Dark2
+
+import pyzebra
+
+
+javaScript = """
+let j = 0;
+for (let i = 0; i < js_data.data['fname'].length; i++) {
+    if (js_data.data['content'][i] === "") continue;
+
+    setTimeout(function() {
+        const blob = new Blob([js_data.data['content'][i]], {type: 'text/plain'})
+        const link = document.createElement('a');
+        document.body.appendChild(link);
+        const url = window.URL.createObjectURL(blob);
+        link.href = url;
+        link.download = js_data.data['fname'][i];
+        link.click();
+        window.URL.revokeObjectURL(url);
+        document.body.removeChild(link);
+    }, 100 * j)
+
+    j++;
+}
+"""
+
+ANG_CHUNK_DEFAULTS = {"2theta": 30, "gamma": 30, "omega": 30, "chi": 35, "phi": 35, "nu": 10}
+SORT_OPT_BI = ["2theta", "chi", "phi", "omega"]
+SORT_OPT_NB = ["gamma", "nu", "omega"]
+
+
+def create():
+    ang_lims = None
+    cif_data = None
+    params = None
+    res_files = {}
+    js_data = ColumnDataSource(data=dict(content=[""], fname=[""]))
+
+    anglim_div = Div(text="Angular min/max limits:", margin=(5, 5, 0, 5))
+    sttgamma_ti = TextInput(title="stt/gamma", width=100)
+    omega_ti = TextInput(title="omega", width=100)
+    chinu_ti = TextInput(title="chi/nu", width=100)
+    phi_ti = TextInput(title="phi", width=100)
+
+    def _update_ang_lims(ang_lims):
+        sttgamma_ti.value = " ".join(ang_lims["gamma"][:2])
+        omega_ti.value = " ".join(ang_lims["omega"][:2])
+        if ang_lims["geom"] == "nb":
+            chinu_ti.value = " ".join(ang_lims["nu"][:2])
+            phi_ti.value = ""
+        else:  # ang_lims["geom"] == "bi"
+            chinu_ti.value = " ".join(ang_lims["chi"][:2])
+            phi_ti.value = " ".join(ang_lims["phi"][:2])
+
+    def _update_params(params):
+        if "WAVE" in params:
+            wavelen_input.value = params["WAVE"]
+        if "SPGR" in params:
+            cryst_space_group.value = params["SPGR"]
+        if "CELL" in params:
+            cryst_cell.value = params["CELL"]
+        if "UBMAT" in params:
+            ub_matrix.value = " ".join(params["UBMAT"])
+        if "HLIM" in params:
+            ranges_hkl.value = params["HLIM"]
+        if "SRANG" in params:
+            ranges_srang.value = params["SRANG"]
+        if "lattiCE" in params:
+            magstruct_lattice.value = params["lattiCE"]
+        if "kvect" in params:
+            magstruct_kvec.value = params["kvect"]
+
+    def open_geom_callback(_attr, _old, new):
+        nonlocal ang_lims
+        with io.StringIO(base64.b64decode(new).decode()) as fileobj:
+            ang_lims = pyzebra.read_geom_file(fileobj)
+        _update_ang_lims(ang_lims)
+
+    open_geom_div = Div(text="Open GEOM:")
+    open_geom = FileInput(accept=".geom", width=200)
+    open_geom.on_change("value", open_geom_callback)
+
+    def open_cfl_callback(_attr, _old, new):
+        nonlocal params
+        with io.StringIO(base64.b64decode(new).decode()) as fileobj:
+            params = pyzebra.read_cfl_file(fileobj)
+            _update_params(params)
+
+    open_cfl_div = Div(text="Open CFL:")
+    open_cfl = FileInput(accept=".cfl", width=200)
+    open_cfl.on_change("value", open_cfl_callback)
+
+    def open_cif_callback(_attr, _old, new):
+        nonlocal cif_data
+        with io.StringIO(base64.b64decode(new).decode()) as fileobj:
+            cif_data = pyzebra.read_cif_file(fileobj)
+            _update_params(cif_data)
+
+    open_cif_div = Div(text="Open CIF:")
+    open_cif = FileInput(accept=".cif", width=200)
+    open_cif.on_change("value", open_cif_callback)
+
+    wavelen_div = Div(text="Wavelength:", margin=(5, 5, 0, 5))
+    wavelen_input = TextInput(title="value", width=70)
+
+    def wavelen_select_callback(_attr, _old, new):
+        if new:
+            wavelen_input.value = new
+        else:
+            wavelen_input.value = ""
+
+    wavelen_select = Select(
+        title="preset", options=["", "0.788", "1.178", "1.383", "2.305"], width=70
+    )
+    wavelen_select.on_change("value", wavelen_select_callback)
+
+    cryst_div = Div(text="Crystal structure:", margin=(5, 5, 0, 5))
+    cryst_space_group = TextInput(title="space group", width=100)
+    cryst_cell = TextInput(title="cell", width=250)
+
+    def ub_matrix_calc_callback():
+        params = dict()
+        params["SPGR"] = cryst_space_group.value
+        params["CELL"] = cryst_cell.value
+        ub = pyzebra.calc_ub_matrix(params)
+        ub_matrix.value = " ".join(ub)
+
+    ub_matrix_calc = Button(label="UB matrix:", button_type="primary", width=100)
+    ub_matrix_calc.on_click(ub_matrix_calc_callback)
+
+    ub_matrix = TextInput(title="\u200B", width=600)
+
+    ranges_div = Div(text="Ranges:", margin=(5, 5, 0, 5))
+    ranges_hkl = TextInput(title="HKL", value="-25 25 -25 25 -25 25", width=250)
+    ranges_srang = TextInput(title="sin(​θ​)/λ", value="0.0 0.7", width=100)
+
+    magstruct_div = Div(text="Magnetic structure:", margin=(5, 5, 0, 5))
+    magstruct_lattice = TextInput(title="lattice", width=100)
+    magstruct_kvec = TextAreaInput(title="k vector", width=150)
+
+    def sorting0_callback(_attr, _old, new):
+        sorting_0_dt.value = ANG_CHUNK_DEFAULTS[new]
+
+    def sorting1_callback(_attr, _old, new):
+        sorting_1_dt.value = ANG_CHUNK_DEFAULTS[new]
+
+    def sorting2_callback(_attr, _old, new):
+        sorting_2_dt.value = ANG_CHUNK_DEFAULTS[new]
+
+    sorting_0 = Select(title="1st", width=100)
+    sorting_0.on_change("value", sorting0_callback)
+    sorting_0_dt = NumericInput(title="Δ", width=70)
+    sorting_1 = Select(title="2nd", width=100)
+    sorting_1.on_change("value", sorting1_callback)
+    sorting_1_dt = NumericInput(title="Δ", width=70)
+    sorting_2 = Select(title="3rd", width=100)
+    sorting_2.on_change("value", sorting2_callback)
+    sorting_2_dt = NumericInput(title="Δ", width=70)
+
+    def geom_radiogroup_callback(_attr, _old, new):
+        nonlocal ang_lims, params
+        if new == 0:
+            geom_file = pyzebra.get_zebraBI_default_geom_file()
+            sort_opt = SORT_OPT_BI
+        else:
+            geom_file = pyzebra.get_zebraNB_default_geom_file()
+            sort_opt = SORT_OPT_NB
+        cfl_file = pyzebra.get_zebra_default_cfl_file()
+
+        ang_lims = pyzebra.read_geom_file(geom_file)
+        _update_ang_lims(ang_lims)
+        params = pyzebra.read_cfl_file(cfl_file)
+        _update_params(params)
+
+        sorting_0.options = sorting_1.options = sorting_2.options = sort_opt
+        sorting_0.value = sort_opt[0]
+        sorting_1.value = sort_opt[1]
+        sorting_2.value = sort_opt[2]
+
+    geom_radiogroup_div = Div(text="Geometry:", margin=(5, 5, 0, 5))
+    geom_radiogroup = RadioGroup(labels=["bisecting", "normal beam"], width=150)
+    geom_radiogroup.on_change("active", geom_radiogroup_callback)
+    geom_radiogroup.active = 0
+
+    def go_button_callback():
+        ang_lims["gamma"][0], ang_lims["gamma"][1] = sttgamma_ti.value.strip().split()
+        ang_lims["omega"][0], ang_lims["omega"][1] = omega_ti.value.strip().split()
+        if ang_lims["geom"] == "nb":
+            ang_lims["nu"][0], ang_lims["nu"][1] = chinu_ti.value.strip().split()
+        else:  # ang_lims["geom"] == "bi"
+            ang_lims["chi"][0], ang_lims["chi"][1] = chinu_ti.value.strip().split()
+            ang_lims["phi"][0], ang_lims["phi"][1] = phi_ti.value.strip().split()
+
+        if cif_data:
+            params.update(cif_data)
+
+        params["WAVE"] = wavelen_input.value
+        params["SPGR"] = cryst_space_group.value
+        params["CELL"] = cryst_cell.value
+        params["UBMAT"] = ub_matrix.value.split()
+        params["HLIM"] = ranges_hkl.value
+        params["SRANG"] = ranges_srang.value
+        params["lattiCE"] = magstruct_lattice.value
+        kvects = magstruct_kvec.value.split("\n")
+
+        with tempfile.TemporaryDirectory() as temp_dir:
+            geom_path = os.path.join(temp_dir, "zebra.geom")
+            if open_geom.value:
+                geom_template = io.StringIO(base64.b64decode(open_geom.value).decode())
+            else:
+                geom_template = None
+            pyzebra.export_geom_file(geom_path, ang_lims, geom_template)
+
+            print(f"Content of {geom_path}:")
+            with open(geom_path) as f:
+                print(f.read())
+
+            priority = [sorting_0.value, sorting_1.value, sorting_2.value]
+            chunks = [sorting_0_dt.value, sorting_1_dt.value, sorting_2_dt.value]
+            if geom_radiogroup.active == 0:
+                sort_hkl_file = pyzebra.sort_hkl_file_bi
+                priority.extend(set(SORT_OPT_BI) - set(priority))
+            else:
+                sort_hkl_file = pyzebra.sort_hkl_file_nb
+
+            # run sxtal_refgen for each kvect provided
+            for i, kvect in enumerate(kvects, start=1):
+                params["kvect"] = kvect
+                if open_cfl.filename:
+                    base_fname = f"{os.path.splitext(open_cfl.filename)[0]}_{i}"
+                else:
+                    base_fname = f"zebra_{i}"
+
+                cfl_path = os.path.join(temp_dir, base_fname + ".cfl")
+                if open_cfl.value:
+                    cfl_template = io.StringIO(base64.b64decode(open_cfl.value).decode())
+                else:
+                    cfl_template = None
+                pyzebra.export_cfl_file(cfl_path, params, cfl_template)
+
+                print(f"Content of {cfl_path}:")
+                with open(cfl_path) as f:
+                    print(f.read())
+
+                comp_proc = subprocess.run(
+                    [pyzebra.SXTAL_REFGEN_PATH, cfl_path],
+                    cwd=temp_dir,
+                    check=True,
+                    stdout=subprocess.PIPE,
+                    stderr=subprocess.STDOUT,
+                    text=True,
+                )
+                print(" ".join(comp_proc.args))
+                print(comp_proc.stdout)
+
+                if i == 1:  # all hkl files are identical, so keep only one
+                    hkl_fname = base_fname + ".hkl"
+                    hkl_fpath = os.path.join(temp_dir, hkl_fname)
+                    with open(hkl_fpath) as f:
+                        res_files[hkl_fname] = f.read()
+
+                    hkl_fname_sorted = base_fname + "_sorted.hkl"
+                    hkl_fpath_sorted = os.path.join(temp_dir, hkl_fname_sorted)
+                    sort_hkl_file(hkl_fpath, hkl_fpath_sorted, priority, chunks)
+                    with open(hkl_fpath_sorted) as f:
+                        res_files[hkl_fname_sorted] = f.read()
+
+                mhkl_fname = base_fname + ".mhkl"
+                mhkl_fpath = os.path.join(temp_dir, mhkl_fname)
+                with open(mhkl_fpath) as f:
+                    res_files[mhkl_fname] = f.read()
+
+                mhkl_fname_sorted = base_fname + "_sorted.mhkl"
+                mhkl_fpath_sorted = os.path.join(temp_dir, hkl_fname_sorted)
+                sort_hkl_file(mhkl_fpath, mhkl_fpath_sorted, priority, chunks)
+                with open(mhkl_fpath_sorted) as f:
+                    res_files[mhkl_fname_sorted] = f.read()
+
+            created_lists.options = list(res_files)
+
+    go_button = Button(label="GO", button_type="primary", width=50)
+    go_button.on_click(go_button_callback)
+
+    def created_lists_callback(_attr, _old, new):
+        sel_file = new[0]
+        file_text = res_files[sel_file]
+        preview_lists.value = file_text
+        js_data.data.update(content=[file_text], fname=[sel_file])
+
+    created_lists = MultiSelect(title="Created lists:", width=200, height=150)
+    created_lists.on_change("value", created_lists_callback)
+    preview_lists = TextAreaInput(title="Preview selected list:", width=600, height=150)
+
+    download_file = Button(label="Download file", button_type="success", width=200)
+    download_file.js_on_click(CustomJS(args={"js_data": js_data}, code=javaScript))
+    plot_list = Button(label="Plot selected list", button_type="primary", width=200, disabled=True)
+
+    measured_data_div = Div(text="Measured data:")
+    measured_data = FileInput(accept=".ccl", multiple=True, width=200)
+
+    min_grid_x = -10
+    max_grid_x = 10
+    min_grid_y = -5
+    max_grid_y = 5
+    cmap = Dark2[8]
+    syms = ["circle", "inverted_triangle", "square", "diamond", "star", "triangle"]
+
+    # Define resolution function
+    def _res_fun(stt, wave, res_mult):
+        expr = np.tan(stt / 2 * np.pi / 180)
+        fwhm = np.sqrt(0.4639 * expr ** 2 - 0.4452 * expr + 0.1506) * res_mult  # res in deg
+        return fwhm
+
+    def plot_file_callback():
+        orth_dir = list(map(float, hkl_normal.value.split()))
+        cut_tol = hkl_delta.value
+        cut_or = hkl_cut.value
+        x_dir = list(map(float, hkl_in_plane_x.value.split()))
+        y_dir = list(map(float, hkl_in_plane_y.value.split()))
+
+        k = np.array(k_vectors.value.split()).astype(float).reshape(3, 3)
+        tol_k = 0.1
+
+        # Plotting options
+        grid_flag = 1
+        grid_minor_flag = 1
+        grid_div = 2  # Number of minor division lines per unit
+
+        # different symbols based on file number
+        file_flag = 0 in disting_opt_cb.active
+        # scale marker size according to intensity
+        intensity_flag = 1 in disting_opt_cb.active
+        # use color to mark different propagation vectors
+        prop_legend_flag = 2 in disting_opt_cb.active
+        # use resolution ellipsis
+        res_flag = disting_opt_rb.active
+        # multiplier for resolution function (in case of samples with large mosaicity)
+        res_mult = res_mult_ni.value
+
+        md_fnames = measured_data.filename
+        md_fdata = measured_data.value
+
+        # Load first data cile, read angles and define matrices to perform conversion to cartesian coordinates and back
+        with io.StringIO(base64.b64decode(md_fdata[0]).decode()) as file:
+            _, ext = os.path.splitext(md_fnames[0])
+            try:
+                file_data = pyzebra.parse_1D(file, ext)
+            except:
+                print(f"Error loading {md_fnames[0]}")
+                return
+
+        alpha = file_data[0]["alpha_cell"] * np.pi / 180.0
+        beta = file_data[0]["beta_cell"] * np.pi / 180.0
+        gamma = file_data[0]["gamma_cell"] * np.pi / 180.0
+
+        # reciprocal angle parameters
+        alpha_star = np.arccos(
+            (np.cos(beta) * np.cos(gamma) - np.cos(alpha)) / (np.sin(beta) * np.sin(gamma))
+        )
+        beta_star = np.arccos(
+            (np.cos(alpha) * np.cos(gamma) - np.cos(beta)) / (np.sin(alpha) * np.sin(gamma))
+        )
+        gamma_star = np.arccos(
+            (np.cos(alpha) * np.cos(beta) - np.cos(gamma)) / (np.sin(alpha) * np.sin(beta))
+        )
+
+        # conversion matrix:
+        M = np.array(
+            [
+                [1, np.cos(gamma_star), np.cos(beta_star)],
+                [0, np.sin(gamma_star), -np.sin(beta_star) * np.cos(alpha)],
+                [0, 0, np.sin(beta_star) * np.sin(alpha)],
+            ]
+        )
+        M_inv = np.linalg.inv(M)
+
+        # Calculate in-plane y-direction
+        x_c = M @ x_dir
+        y_c = M @ y_dir
+        o_c = M @ orth_dir
+
+        # Normalize all directions
+        y_c = y_c / np.linalg.norm(y_c)
+        x_c = x_c / np.linalg.norm(x_c)
+        o_c = o_c / np.linalg.norm(o_c)
+
+        # Read all data
+        hkl_coord = []
+        intensity_vec = []
+        k_flag_vec = []
+        file_flag_vec = []
+        res_vec_x = []
+        res_vec_y = []
+        res_N = 10
+
+        for j in range(len(md_fnames)):
+            with io.StringIO(base64.b64decode(md_fdata[j]).decode()) as file:
+                _, ext = os.path.splitext(md_fnames[j])
+                try:
+                    file_data = pyzebra.parse_1D(file, ext)
+                except:
+                    print(f"Error loading {md_fnames[j]}")
+                    return
+
+            # Loop throguh all data
+            for scan in file_data:
+                om = scan["omega"]
+                gammad = scan["twotheta"]
+                chi = scan["chi"]
+                phi = scan["phi"]
+                nud = 0  # 1d detector
+                ub = scan["ub"]
+                ddist = float(scan["detectorDistance"])
+                counts = scan["counts"]
+                mon = scan["monitor"]
+
+                # Determine wavelength from mcvl value (is wavelength stored anywhere???)
+                mcvl = scan["mcvl"]
+                if mcvl == 2.2:
+                    wave = 1.178
+                elif mcvl == 7.0:
+                    wave = 1.383
+                else:
+                    wave = 2.3
+
+                # Calculate resolution in degrees
+                res = _res_fun(gammad, wave, res_mult)
+
+                # convert to resolution in hkl along scan line
+                ang2hkl_1d = pyzebra.ang2hkl_1d
+                res_x = []
+                res_y = []
+                for _om in np.linspace(om[0], om[-1], num=res_N):
+                    expr1 = ang2hkl_1d(wave, ddist, gammad, _om + res / 2, chi, phi, nud, ub)
+                    expr2 = ang2hkl_1d(wave, ddist, gammad, _om - res / 2, chi, phi, nud, ub)
+                    hkl_temp = M @ (np.abs(expr1 - expr2) / 2)
+                    res_x.append(hkl_temp[0])
+                    res_y.append(hkl_temp[1])
+
+                # Get first and final hkl
+                hkl1 = ang2hkl_1d(wave, ddist, gammad, om[0], chi, phi, nud, ub)
+                hkl2 = ang2hkl_1d(wave, ddist, gammad, om[-1], chi, phi, nud, ub)
+
+                # Get hkl at best intensity
+                hkl_m = ang2hkl_1d(wave, ddist, gammad, om[np.argmax(counts)], chi, phi, nud, ub)
+
+                # Estimate intensity for marker size scaling
+                y1 = counts[0]
+                y2 = counts[-1]
+                x1 = om[0]
+                x2 = om[-1]
+                a = (y1 - y2) / (x1 - x2)
+                b = y1 - a * x1
+                intensity_exp = np.sum(counts - (a * om + b))
+                c = int(intensity_exp / mon * 10000)
+
+                # Recognize k_flag_vec
+                min_hkl_m = np.minimum(1 - hkl_m % 1, hkl_m % 1)
+                for j2, _k in enumerate(k):
+                    if all(np.abs(min_hkl_m - _k) < tol_k):
+                        k_flag_vec.append(j2)
+                        break
+                else:
+                    k_flag_vec.append(len(k))
+
+                # Save data
+                hkl_coord.append([hkl1, hkl2, hkl_m])
+                intensity_vec.append(c)
+                file_flag_vec.append(j)
+                res_vec_x.append(res_x)
+                res_vec_y.append(res_y)
+
+        plot.x_range.start = plot.x_range.reset_start = -2
+        plot.x_range.end = plot.x_range.reset_end = 5
+        plot.y_range.start = plot.y_range.reset_start = -4
+        plot.y_range.end = plot.y_range.reset_end = 3.5
+
+        xs, ys = [], []
+        xs_minor, ys_minor = [], []
+        if grid_flag:
+            for yy in np.arange(min_grid_y, max_grid_y, 1):
+                hkl1 = M @ [0, yy, 0]
+                xs.append([min_grid_y, max_grid_y])
+                ys.append([hkl1[1], hkl1[1]])
+
+            for xx in np.arange(min_grid_x, max_grid_x, 1):
+                hkl1 = M @ [xx, min_grid_x, 0]
+                hkl2 = M @ [xx, max_grid_x, 0]
+                xs.append([hkl1[0], hkl2[0]])
+                ys.append([hkl1[1], hkl2[1]])
+
+            if grid_minor_flag:
+                for yy in np.arange(min_grid_y, max_grid_y, 1 / grid_div):
+                    hkl1 = M @ [0, yy, 0]
+                    xs_minor.append([min_grid_y, max_grid_y])
+                    ys_minor.append([hkl1[1], hkl1[1]])
+
+                for xx in np.arange(min_grid_x, max_grid_x, 1 / grid_div):
+                    hkl1 = M @ [xx, min_grid_x, 0]
+                    hkl2 = M @ [xx, max_grid_x, 0]
+                    xs_minor.append([hkl1[0], hkl2[0]])
+                    ys_minor.append([hkl1[1], hkl2[1]])
+
+        grid_source.data.update(xs=xs, ys=ys)
+        minor_grid_source.data.update(xs=xs_minor, ys=ys_minor)
+
+        el_x, el_y, el_w, el_h, el_c = [], [], [], [], []
+        scan_xs, scan_ys, scan_x, scan_y = [], [], [], []
+        scan_m, scan_s, scan_c, scan_l = [], [], [], []
+        for j in range(len(hkl_coord)):
+            # Get middle hkl from list
+            hklm = M @ hkl_coord[j][2]
+
+            # Decide if point is in the cut
+            proj = np.dot(hklm, o_c)
+            if abs(proj - cut_or) >= cut_tol:
+                continue
+
+            hkl1 = M @ hkl_coord[j][0]
+            hkl2 = M @ hkl_coord[j][1]
+
+            if intensity_flag:
+                markersize = max(1, int(intensity_vec[j] / max(intensity_vec) * 20))
+            else:
+                markersize = 4
+
+            if file_flag:
+                plot_symbol = syms[file_flag_vec[j]]
+            else:
+                plot_symbol = "circle"
+
+            if prop_legend_flag:
+                col_value = cmap[k_flag_vec[j]]
+            else:
+                col_value = "black"
+
+            if res_flag:
+                # Generate series of ellipses along scan line
+                el_x.extend(np.linspace(hkl1[0], hkl2[0], num=res_N))
+                el_y.extend(np.linspace(hkl1[1], hkl2[1], num=res_N))
+                el_w.extend(np.array(res_vec_x[j]) * 2)
+                el_h.extend(np.array(res_vec_y[j]) * 2)
+                el_c.extend([col_value] * res_N)
+            else:
+                # Plot scan line
+                scan_xs.append([hkl1[0], hkl2[0]])
+                scan_ys.append([hkl1[1], hkl2[1]])
+
+                # Plot middle point of scan
+                scan_x.append(hklm[0])
+                scan_y.append(hklm[1])
+                scan_m.append(plot_symbol)
+                scan_s.append(markersize)
+
+                # Color and legend label
+                scan_c.append(col_value)
+                scan_l.append(md_fnames[file_flag_vec[j]])
+
+        ellipse_source.data.update(x=el_x, y=el_y, w=el_w, h=el_h, c=el_c)
+        scan_source.data.update(
+            xs=scan_xs, ys=scan_ys, x=scan_x, y=scan_y, m=scan_m, s=scan_s, c=scan_c, l=scan_l,
+        )
+
+        arrow1.visible = True
+        arrow1.x_end = x_c[0]
+        arrow1.y_end = x_c[1]
+        arrow2.visible = True
+        arrow2.x_end = y_c[0]
+        arrow2.y_end = y_c[1]
+
+        kvect_source.data.update(
+            text_x=[x_c[0] / 2, y_c[0] / 2 - 0.1],
+            text_y=[x_c[1] - 0.1, y_c[1] / 2],
+            text=["h", "k"],
+        )
+
+        # Legend items for different file entries (symbol)
+        legend_items = []
+        if not res_flag and file_flag:
+            labels, inds = np.unique(scan_source.data["l"], return_index=True)
+            for label, ind in zip(labels, inds):
+                legend_items.append(LegendItem(label=label, renderers=[scatter], index=ind))
+
+        # Legend items for propagation vector (color)
+        if prop_legend_flag:
+            if res_flag:
+                source, render = ellipse_source, ellipse
+            else:
+                source, render = scan_source, mline
+
+            labels, inds = np.unique(source.data["c"], return_index=True)
+            for label, ind in zip(labels, inds):
+                label = f"k={k[cmap.index(label)]}"
+                legend_items.append(LegendItem(label=label, renderers=[render], index=ind))
+
+        plot.legend.items = legend_items
+
+    plot_file = Button(label="Plot selected file(s)", button_type="primary", width=200)
+    plot_file.on_click(plot_file_callback)
+
+    plot = Plot(x_range=Range1d(), y_range=Range1d(), plot_height=450, plot_width=600)
+    plot.add_tools(PanTool(), WheelZoomTool(), BoxZoomTool(), ResetTool())
+    plot.toolbar.logo = None
+
+    plot.add_layout(LinearAxis(), place="left")
+    plot.add_layout(LinearAxis(), place="below")
+
+    arrow1 = Arrow(x_start=0, y_start=0, x_end=0, y_end=0, end=NormalHead(size=10), visible=False)
+    plot.add_layout(arrow1)
+    arrow2 = Arrow(x_start=0, y_start=0, x_end=0, y_end=0, end=NormalHead(size=10), visible=False)
+    plot.add_layout(arrow2)
+
+    kvect_source = ColumnDataSource(dict(text_x=[], text_y=[], text=[]))
+    plot.add_glyph(kvect_source, Text(x="text_x", y="text_y", text="text"))
+
+    grid_source = ColumnDataSource(dict(xs=[], ys=[]))
+    minor_grid_source = ColumnDataSource(dict(xs=[], ys=[]))
+    plot.add_glyph(grid_source, MultiLine(xs="xs", ys="ys", line_color="gray"))
+    plot.add_glyph(
+        minor_grid_source, MultiLine(xs="xs", ys="ys", line_color="gray", line_dash="dotted")
+    )
+
+    ellipse_source = ColumnDataSource(dict(x=[], y=[], w=[], h=[], c=[]))
+    ellipse = plot.add_glyph(
+        ellipse_source, Ellipse(x="x", y="y", width="w", height="h", fill_color="c", line_color="c")
+    )
+
+    scan_source = ColumnDataSource(dict(xs=[], ys=[], x=[], y=[], m=[], s=[], c=[], l=[]))
+    mline = plot.add_glyph(scan_source, MultiLine(xs="xs", ys="ys", line_color="c"))
+    scatter = plot.add_glyph(
+        scan_source, Scatter(x="x", y="y", marker="m", size="s", fill_color="c", line_color="c")
+    )
+
+    plot.add_layout(Legend(items=[], location="top_left", click_policy="hide"))
+
+    hkl_div = Div(text="HKL:", margin=(5, 5, 0, 5))
+    hkl_normal = TextInput(title="normal", value="0 0 1", width=70)
+    hkl_cut = Spinner(title="cut", value=0, step=0.1, width=70)
+    hkl_delta = NumericInput(title="delta", value=0.1, mode="float", width=70)
+    hkl_in_plane_x = TextInput(title="in-plane X", value="1 0 0", width=70)
+    hkl_in_plane_y = TextInput(title="in-plane Y", value="0 1 0", width=70)
+
+    disting_opt_div = Div(text="Distinguish options:", margin=(5, 5, 0, 5))
+    disting_opt_cb = CheckboxGroup(
+        labels=["files (symbols)", "intensities (size)", "k vectors nucl/magn (colors)"],
+        active=[0, 1, 2],
+        width=200,
+    )
+    disting_opt_rb = RadioGroup(
+        labels=["scan direction", "resolution ellipsoid"], active=0, width=200
+    )
+
+    k_vectors = TextAreaInput(
+        title="k vectors:", value="0.0 0.0 0.0\n0.5 0.0 0.0\n0.5 0.5 0.0", width=150,
+    )
+    res_mult_ni = NumericInput(title="Resolution mult:", value=10, mode="int", width=100)
+
+    fileinput_layout = row(open_cfl_div, open_cfl, open_cif_div, open_cif, open_geom_div, open_geom)
+
+    geom_layout = column(geom_radiogroup_div, geom_radiogroup)
+    wavelen_layout = column(wavelen_div, row(wavelen_select, wavelen_input))
+    anglim_layout = column(anglim_div, row(sttgamma_ti, omega_ti, chinu_ti, phi_ti))
+    cryst_layout = column(cryst_div, row(cryst_space_group, cryst_cell))
+    ubmat_layout = row(column(Spacer(height=18), ub_matrix_calc), ub_matrix)
+    ranges_layout = column(ranges_div, row(ranges_hkl, ranges_srang))
+    magstruct_layout = column(magstruct_div, row(magstruct_lattice, magstruct_kvec))
+    sorting_layout = row(
+        sorting_0,
+        sorting_0_dt,
+        Spacer(width=30),
+        sorting_1,
+        sorting_1_dt,
+        Spacer(width=30),
+        sorting_2,
+        sorting_2_dt,
+    )
+
+    column1_layout = column(
+        fileinput_layout,
+        Spacer(height=10),
+        row(geom_layout, wavelen_layout, Spacer(width=50), anglim_layout),
+        cryst_layout,
+        ubmat_layout,
+        row(ranges_layout, Spacer(width=50), magstruct_layout),
+        row(sorting_layout, Spacer(width=30), column(Spacer(height=18), go_button)),
+        row(created_lists, preview_lists),
+        row(download_file, plot_list),
+    )
+
+    hkl_layout = column(
+        hkl_div,
+        row(hkl_normal, hkl_cut, hkl_delta, Spacer(width=10), hkl_in_plane_x, hkl_in_plane_y),
+    )
+    disting_layout = column(disting_opt_div, row(disting_opt_cb, disting_opt_rb))
+
+    column2_layout = column(
+        row(measured_data_div, measured_data, plot_file),
+        plot,
+        row(hkl_layout, k_vectors),
+        row(disting_layout, res_mult_ni),
+    )
+
+    tab_layout = row(column1_layout, column2_layout)
+
+    return Panel(child=tab_layout, title="ccl prepare")

pyzebra/app/panel_hdf_anatric.py

@@ -1,5 +1,6 @@
 import base64
 import io
+import os
 import re
 import tempfile
 
@@ -10,8 +11,9 @@ from bokeh.models import (
     Div,
     FileInput,
     Panel,
-    RadioButtonGroup,
     Select,
+    Spacer,
+    Tabs,
     TextAreaInput,
     TextInput,
 )
@@ -28,7 +30,7 @@ def create():
         config.load_from_file(file)
 
         logfile_textinput.value = config.logfile
-        logfile_verbosity_select.value = config.logfile_verbosity
+        logfile_verbosity.value = config.logfile_verbosity
 
         filelist_type.value = config.filelist_type
         filelist_format_textinput.value = config.filelist_format
@@ -43,11 +45,16 @@ def create():
         ub_textareainput.value = config.crystal_UB
 
         dataFactory_implementation_select.value = config.dataFactory_implementation
-        dataFactory_dist1_textinput.value = config.dataFactory_dist1
+        if config.dataFactory_dist1 is not None:
+            dataFactory_dist1_textinput.value = config.dataFactory_dist1
+        if config.dataFactory_dist2 is not None:
+            dataFactory_dist2_textinput.value = config.dataFactory_dist2
+        if config.dataFactory_dist3 is not None:
+            dataFactory_dist3_textinput.value = config.dataFactory_dist3
         reflectionPrinter_format_select.value = config.reflectionPrinter_format
 
-        set_active_widgets(config.algorithm)
         if config.algorithm == "adaptivemaxcog":
+            algorithm_params.active = 0
             threshold_textinput.value = config.threshold
             shell_textinput.value = config.shell
             steepness_textinput.value = config.steepness
@@ -56,6 +63,7 @@ def create():
             aps_window_textinput.value = str(tuple(map(int, config.aps_window.values())))
 
         elif config.algorithm == "adaptivedynamic":
+            algorithm_params.active = 1
             adm_window_textinput.value = str(tuple(map(int, config.adm_window.values())))
             border_textinput.value = str(tuple(map(int, config.border.values())))
             minWindow_textinput.value = str(tuple(map(int, config.minWindow.values())))
@@ -65,46 +73,16 @@ def create():
             loop_textinput.value = config.loop
             minPeakCount_textinput.value = config.minPeakCount
             displacementCurve_textinput.value = "\n".join(map(str, config.displacementCurve))
+
         else:
             raise ValueError("Unknown processing mode.")
 
-    def set_active_widgets(implementation):
-        if implementation == "adaptivemaxcog":
-            mode_radio_button_group.active = 0
-            disable_adaptivemaxcog = False
-            disable_adaptivedynamic = True
-
-        elif implementation == "adaptivedynamic":
-            mode_radio_button_group.active = 1
-            disable_adaptivemaxcog = True
-            disable_adaptivedynamic = False
-        else:
-            raise ValueError("Implementation can be either 'adaptivemaxcog' or 'adaptivedynamic'")
-
-        threshold_textinput.disabled = disable_adaptivemaxcog
-        shell_textinput.disabled = disable_adaptivemaxcog
-        steepness_textinput.disabled = disable_adaptivemaxcog
-        duplicateDistance_textinput.disabled = disable_adaptivemaxcog
-        maxequal_textinput.disabled = disable_adaptivemaxcog
-        aps_window_textinput.disabled = disable_adaptivemaxcog
-
-        adm_window_textinput.disabled = disable_adaptivedynamic
-        border_textinput.disabled = disable_adaptivedynamic
-        minWindow_textinput.disabled = disable_adaptivedynamic
-        reflectionFile_textinput.disabled = disable_adaptivedynamic
-        targetMonitor_textinput.disabled = disable_adaptivedynamic
-        smoothSize_textinput.disabled = disable_adaptivedynamic
-        loop_textinput.disabled = disable_adaptivedynamic
-        minPeakCount_textinput.disabled = disable_adaptivedynamic
-        displacementCurve_textinput.disabled = disable_adaptivedynamic
-
-    upload_div = Div(text="Open XML configuration file:")
-
     def upload_button_callback(_attr, _old, new):
         with io.BytesIO(base64.b64decode(new)) as file:
             _load_config_file(file)
 
-    upload_button = FileInput(accept=".xml")
+    upload_div = Div(text="Open .xml config:")
+    upload_button = FileInput(accept=".xml", width=200)
     upload_button.on_change("value", upload_button_callback)
 
     # General parameters
@@ -112,16 +90,14 @@ def create():
     def logfile_textinput_callback(_attr, _old, new):
         config.logfile = new
 
-    logfile_textinput = TextInput(title="Logfile:", value="logfile.log", width=520)
+    logfile_textinput = TextInput(title="Logfile:", value="logfile.log")
     logfile_textinput.on_change("value", logfile_textinput_callback)
 
-    def logfile_verbosity_select_callback(_attr, _old, new):
+    def logfile_verbosity_callback(_attr, _old, new):
         config.logfile_verbosity = new
 
-    logfile_verbosity_select = Select(
-        title="verbosity:", options=["0", "5", "10", "15", "30"], width=70
-    )
-    logfile_verbosity_select.on_change("value", logfile_verbosity_select_callback)
+    logfile_verbosity = TextInput(title="verbosity:", width=70)
+    logfile_verbosity.on_change("value", logfile_verbosity_callback)
 
     # ---- FileList
     def filelist_type_callback(_attr, _old, new):
@@ -133,7 +109,7 @@ def create():
     def filelist_format_textinput_callback(_attr, _old, new):
         config.filelist_format = new
 
-    filelist_format_textinput = TextInput(title="format:", width=490)
+    filelist_format_textinput = TextInput(title="format:", width=290)
     filelist_format_textinput.on_change("value", filelist_format_textinput_callback)
 
     def filelist_datapath_textinput_callback(_attr, _old, new):
@@ -148,20 +124,20 @@ def create():
             ranges.append(re.findall(r"\b\d+\b", line))
         config.filelist_ranges = ranges
 
-    filelist_ranges_textareainput = TextAreaInput(title="ranges:", height=100)
+    filelist_ranges_textareainput = TextAreaInput(title="ranges:", rows=1)
     filelist_ranges_textareainput.on_change("value", filelist_ranges_textareainput_callback)
 
     # ---- crystal
     def crystal_sample_textinput_callback(_attr, _old, new):
         config.crystal_sample = new
 
-    crystal_sample_textinput = TextInput(title="Sample Name:")
+    crystal_sample_textinput = TextInput(title="Sample Name:", width=290)
     crystal_sample_textinput.on_change("value", crystal_sample_textinput_callback)
 
     def lambda_textinput_callback(_attr, _old, new):
         config.crystal_lambda = new
 
-    lambda_textinput = TextInput(title="lambda:", width=140)
+    lambda_textinput = TextInput(title="lambda:", width=100)
     lambda_textinput.on_change("value", lambda_textinput_callback)
 
     def ub_textareainput_callback(_attr, _old, new):
@@ -173,19 +149,19 @@ def create():
     def zeroOM_textinput_callback(_attr, _old, new):
         config.crystal_zeroOM = new
 
-    zeroOM_textinput = TextInput(title="zeroOM:", width=140)
+    zeroOM_textinput = TextInput(title="zeroOM:", width=100)
     zeroOM_textinput.on_change("value", zeroOM_textinput_callback)
 
     def zeroSTT_textinput_callback(_attr, _old, new):
         config.crystal_zeroSTT = new
 
-    zeroSTT_textinput = TextInput(title="zeroSTT:", width=140)
+    zeroSTT_textinput = TextInput(title="zeroSTT:", width=100)
     zeroSTT_textinput.on_change("value", zeroSTT_textinput_callback)
 
     def zeroCHI_textinput_callback(_attr, _old, new):
         config.crystal_zeroCHI = new
 
-    zeroCHI_textinput = TextInput(title="zeroCHI:", width=140)
+    zeroCHI_textinput = TextInput(title="zeroCHI:", width=100)
     zeroCHI_textinput.on_change("value", zeroCHI_textinput_callback)
 
     # ---- DataFactory
@@ -193,16 +169,28 @@ def create():
         config.dataFactory_implementation = new
 
     dataFactory_implementation_select = Select(
-        title="DataFactory implementation:", options=DATA_FACTORY_IMPLEMENTATION, width=300,
+        title="DataFactory implement.:", options=DATA_FACTORY_IMPLEMENTATION, width=145,
     )
     dataFactory_implementation_select.on_change("value", dataFactory_implementation_select_callback)
 
     def dataFactory_dist1_textinput_callback(_attr, _old, new):
         config.dataFactory_dist1 = new
 
-    dataFactory_dist1_textinput = TextInput(title="dist1:", width=290)
+    dataFactory_dist1_textinput = TextInput(title="dist1:", width=75)
     dataFactory_dist1_textinput.on_change("value", dataFactory_dist1_textinput_callback)
 
+    def dataFactory_dist2_textinput_callback(_attr, _old, new):
+        config.dataFactory_dist2 = new
+
+    dataFactory_dist2_textinput = TextInput(title="dist2:", width=75)
+    dataFactory_dist2_textinput.on_change("value", dataFactory_dist2_textinput_callback)
+
+    def dataFactory_dist3_textinput_callback(_attr, _old, new):
+        config.dataFactory_dist3 = new
+
+    dataFactory_dist3_textinput = TextInput(title="dist3:", width=75)
+    dataFactory_dist3_textinput.on_change("value", dataFactory_dist3_textinput_callback)
+
     # ---- BackgroundProcessor
 
     # ---- DetectorEfficency
@@ -212,7 +200,7 @@ def create():
         config.reflectionPrinter_format = new
 
     reflectionPrinter_format_select = Select(
-        title="ReflectionPrinter format:", options=REFLECTION_PRINTER_FORMATS, width=300,
+        title="ReflectionPrinter format:", options=REFLECTION_PRINTER_FORMATS, width=145,
     )
     reflectionPrinter_format_select.on_change("value", reflectionPrinter_format_select_callback)
 
@@ -221,42 +209,42 @@ def create():
     def threshold_textinput_callback(_attr, _old, new):
         config.threshold = new
 
-    threshold_textinput = TextInput(title="Threshold:")
+    threshold_textinput = TextInput(title="Threshold:", width=145)
     threshold_textinput.on_change("value", threshold_textinput_callback)
 
     # ---- shell
     def shell_textinput_callback(_attr, _old, new):
         config.shell = new
 
-    shell_textinput = TextInput(title="Shell:")
+    shell_textinput = TextInput(title="Shell:", width=145)
     shell_textinput.on_change("value", shell_textinput_callback)
 
     # ---- steepness
     def steepness_textinput_callback(_attr, _old, new):
         config.steepness = new
 
-    steepness_textinput = TextInput(title="Steepness:")
+    steepness_textinput = TextInput(title="Steepness:", width=145)
     steepness_textinput.on_change("value", steepness_textinput_callback)
 
     # ---- duplicateDistance
     def duplicateDistance_textinput_callback(_attr, _old, new):
         config.duplicateDistance = new
 
-    duplicateDistance_textinput = TextInput(title="Duplicate Distance:")
+    duplicateDistance_textinput = TextInput(title="Duplicate Distance:", width=145)
     duplicateDistance_textinput.on_change("value", duplicateDistance_textinput_callback)
 
     # ---- maxequal
     def maxequal_textinput_callback(_attr, _old, new):
         config.maxequal = new
 
-    maxequal_textinput = TextInput(title="Max Equal:")
+    maxequal_textinput = TextInput(title="Max Equal:", width=145)
     maxequal_textinput.on_change("value", maxequal_textinput_callback)
 
     # ---- window
     def aps_window_textinput_callback(_attr, _old, new):
         config.aps_window = dict(zip(("x", "y", "z"), re.findall(r"\b\d+\b", new)))
 
-    aps_window_textinput = TextInput(title="Window (x, y, z):")
+    aps_window_textinput = TextInput(title="Window (x, y, z):", width=145)
     aps_window_textinput.on_change("value", aps_window_textinput_callback)
 
     # Adaptive Dynamic Mask Integration (adaptivedynamic)
@@ -264,56 +252,56 @@ def create():
     def adm_window_textinput_callback(_attr, _old, new):
         config.adm_window = dict(zip(("x", "y", "z"), re.findall(r"\b\d+\b", new)))
 
-    adm_window_textinput = TextInput(title="Window (x, y, z):")
+    adm_window_textinput = TextInput(title="Window (x, y, z):", width=145)
     adm_window_textinput.on_change("value", adm_window_textinput_callback)
 
     # ---- border
     def border_textinput_callback(_attr, _old, new):
         config.border = dict(zip(("x", "y", "z"), re.findall(r"\b\d+\b", new)))
 
-    border_textinput = TextInput(title="Border (x, y, z):")
+    border_textinput = TextInput(title="Border (x, y, z):", width=145)
     border_textinput.on_change("value", border_textinput_callback)
 
     # ---- minWindow
     def minWindow_textinput_callback(_attr, _old, new):
         config.minWindow = dict(zip(("x", "y", "z"), re.findall(r"\b\d+\b", new)))
 
-    minWindow_textinput = TextInput(title="Min Window (x, y, z):")
+    minWindow_textinput = TextInput(title="Min Window (x, y, z):", width=145)
     minWindow_textinput.on_change("value", minWindow_textinput_callback)
 
     # ---- reflectionFile
     def reflectionFile_textinput_callback(_attr, _old, new):
         config.reflectionFile = new
 
-    reflectionFile_textinput = TextInput(title="Reflection File:")
+    reflectionFile_textinput = TextInput(title="Reflection File:", width=145)
     reflectionFile_textinput.on_change("value", reflectionFile_textinput_callback)
 
     # ---- targetMonitor
     def targetMonitor_textinput_callback(_attr, _old, new):
         config.targetMonitor = new
 
-    targetMonitor_textinput = TextInput(title="Target Monitor:")
+    targetMonitor_textinput = TextInput(title="Target Monitor:", width=145)
     targetMonitor_textinput.on_change("value", targetMonitor_textinput_callback)
 
     # ---- smoothSize
     def smoothSize_textinput_callback(_attr, _old, new):
         config.smoothSize = new
 
-    smoothSize_textinput = TextInput(title="Smooth Size:")
+    smoothSize_textinput = TextInput(title="Smooth Size:", width=145)
     smoothSize_textinput.on_change("value", smoothSize_textinput_callback)
 
     # ---- loop
     def loop_textinput_callback(_attr, _old, new):
         config.loop = new
 
-    loop_textinput = TextInput(title="Loop:")
+    loop_textinput = TextInput(title="Loop:", width=145)
     loop_textinput.on_change("value", loop_textinput_callback)
 
     # ---- minPeakCount
     def minPeakCount_textinput_callback(_attr, _old, new):
         config.minPeakCount = new
 
-    minPeakCount_textinput = TextInput(title="Min Peak Count:")
+    minPeakCount_textinput = TextInput(title="Min Peak Count:", width=145)
     minPeakCount_textinput.on_change("value", minPeakCount_textinput_callback)
 
     # ---- displacementCurve
@@ -324,86 +312,82 @@ def create():
         config.displacementCurve = maps
 
     displacementCurve_textinput = TextAreaInput(
-        title="Displacement Curve (twotheta, x, y):", height=100
+        title="Displ. Curve (2θ, x, y):", width=145, height=100
     )
     displacementCurve_textinput.on_change("value", displacementCurve_textinput_callback)
 
-    def mode_radio_button_group_callback(active):
-        if active == 0:
+    def algorithm_tabs_callback(_attr, _old, new):
+        if new == 0:
             config.algorithm = "adaptivemaxcog"
-            set_active_widgets("adaptivemaxcog")
         else:
             config.algorithm = "adaptivedynamic"
-            set_active_widgets("adaptivedynamic")
 
-    mode_radio_button_group = RadioButtonGroup(
-        labels=["Adaptive Peak Detection", "Adaptive Dynamic Integration"], active=0
+    algorithm_params = Tabs(
+        tabs=[
+            Panel(
+                child=column(
+                    row(threshold_textinput, shell_textinput, steepness_textinput),
+                    row(duplicateDistance_textinput, maxequal_textinput, aps_window_textinput),
+                ),
+                title="Peak Search",
+            ),
+            Panel(
+                child=column(
+                    row(adm_window_textinput, border_textinput, minWindow_textinput),
+                    row(reflectionFile_textinput, targetMonitor_textinput, smoothSize_textinput),
+                    row(loop_textinput, minPeakCount_textinput, displacementCurve_textinput),
+                ),
+                title="Dynamic Integration",
+            ),
+        ]
     )
-    mode_radio_button_group.on_click(mode_radio_button_group_callback)
-    set_active_widgets("adaptivemaxcog")
+    algorithm_params.on_change("active", algorithm_tabs_callback)
 
     def process_button_callback():
         with tempfile.TemporaryDirectory() as temp_dir:
-            temp_file = temp_dir + "/temp.xml"
+            temp_file = temp_dir + "/config.xml"
             config.save_as(temp_file)
-            pyzebra.anatric(temp_file, anatric_path=doc.anatric_path)
+            pyzebra.anatric(temp_file, anatric_path=doc.anatric_path, cwd=temp_dir)
 
-            with open(config.logfile) as f_log:
+            with open(os.path.join(temp_dir, config.logfile)) as f_log:
                 output_log.value = f_log.read()
 
+            with open(os.path.join(temp_dir, config.reflectionPrinter_file)) as f_res:
+                output_res.value = f_res.read()
+
     process_button = Button(label="Process", button_type="primary")
     process_button.on_click(process_button_callback)
 
-    output_log = TextAreaInput(title="Logfile output:", height=700, disabled=True)
-    output_config = TextAreaInput(title="Current config:", height=700, width=400, disabled=True)
+    output_log = TextAreaInput(title="Logfile output:", height=320, width=465, disabled=True)
+    output_res = TextAreaInput(title="Result output:", height=320, width=465, disabled=True)
+    output_config = TextAreaInput(title="Current config:", height=320, width=465, disabled=True)
+
+    general_params_layout = column(
+        row(column(Spacer(height=2), upload_div), upload_button),
+        row(logfile_textinput, logfile_verbosity),
+        row(filelist_type, filelist_format_textinput),
+        filelist_datapath_textinput,
+        filelist_ranges_textareainput,
+        row(crystal_sample_textinput, lambda_textinput),
+        ub_textareainput,
+        row(zeroOM_textinput, zeroSTT_textinput, zeroCHI_textinput),
+        row(
+            dataFactory_implementation_select,
+            dataFactory_dist1_textinput,
+            dataFactory_dist2_textinput,
+            dataFactory_dist3_textinput,
+        ),
+        row(reflectionPrinter_format_select),
+    )
 
     tab_layout = row(
-        column(
-            upload_div,
-            upload_button,
-            row(logfile_textinput, logfile_verbosity_select),
-            row(filelist_type, filelist_format_textinput),
-            filelist_datapath_textinput,
-            filelist_ranges_textareainput,
-            crystal_sample_textinput,
-            row(lambda_textinput, zeroOM_textinput, zeroSTT_textinput, zeroCHI_textinput),
-            ub_textareainput,
-            row(dataFactory_implementation_select, dataFactory_dist1_textinput),
-            reflectionPrinter_format_select,
-            process_button,
-        ),
-        column(
-            mode_radio_button_group,
-            row(
-                column(
-                    threshold_textinput,
-                    shell_textinput,
-                    steepness_textinput,
-                    duplicateDistance_textinput,
-                    maxequal_textinput,
-                    aps_window_textinput,
-                ),
-                column(
-                    adm_window_textinput,
-                    border_textinput,
-                    minWindow_textinput,
-                    reflectionFile_textinput,
-                    targetMonitor_textinput,
-                    smoothSize_textinput,
-                    loop_textinput,
-                    minPeakCount_textinput,
-                    displacementCurve_textinput,
-                ),
-            ),
-        ),
-        output_config,
-        output_log,
+        general_params_layout,
+        column(output_config, algorithm_params, row(process_button)),
+        column(output_log, output_res),
     )
 
     async def update_config():
-        config.save_as("debug.xml")
-        with open("debug.xml") as f_config:
-            output_config.value = f_config.read()
+        output_config.value = config.tostring()
 
     doc.add_periodic_callback(update_config, 1000)
 

pyzebra/app/panel_hdf_param_study.py

@@ -0,0 +1,591 @@
+import base64
+import io
+import os
+
+import numpy as np
+from bokeh.io import curdoc
+from bokeh.layouts import column, gridplot, row
+from bokeh.models import (
+    BasicTicker,
+    BoxZoomTool,
+    Button,
+    CellEditor,
+    CheckboxGroup,
+    ColumnDataSource,
+    DataRange1d,
+    DataTable,
+    Div,
+    FileInput,
+    Grid,
+    MultiSelect,
+    NumberEditor,
+    NumberFormatter,
+    Image,
+    LinearAxis,
+    LinearColorMapper,
+    Panel,
+    PanTool,
+    Plot,
+    Range1d,
+    ResetTool,
+    Scatter,
+    Select,
+    Spinner,
+    TableColumn,
+    Tabs,
+    Title,
+    WheelZoomTool,
+)
+from bokeh.palettes import Cividis256, Greys256, Plasma256  # pylint: disable=E0611
+
+import pyzebra
+
+IMAGE_W = 256
+IMAGE_H = 128
+IMAGE_PLOT_W = int(IMAGE_W * 2) + 52
+IMAGE_PLOT_H = int(IMAGE_H * 2) + 27
+
+
+def create():
+    doc = curdoc()
+    dataset = []
+    cami_meta = {}
+
+    num_formatter = NumberFormatter(format="0.00", nan_format="")
+
+    def file_select_update():
+        if data_source.value == "proposal number":
+            proposal_path = proposal_textinput.name
+            if proposal_path:
+                file_list = []
+                for file in os.listdir(proposal_path):
+                    if file.endswith(".hdf"):
+                        file_list.append((os.path.join(proposal_path, file), file))
+                file_select.options = file_list
+            else:
+                file_select.options = []
+
+        else:  # "cami file"
+            if not cami_meta:
+                file_select.options = []
+                return
+
+            file_list = cami_meta["filelist"]
+            file_select.options = [(entry, os.path.basename(entry)) for entry in file_list]
+
+    def data_source_callback(_attr, _old, _new):
+        file_select_update()
+
+    data_source = Select(
+        title="Data Source:",
+        value="proposal number",
+        options=["proposal number", "cami file"],
+        width=210,
+    )
+    data_source.on_change("value", data_source_callback)
+
+    doc.add_periodic_callback(file_select_update, 5000)
+
+    def proposal_textinput_callback(_attr, _old, _new):
+        file_select_update()
+
+    proposal_textinput = doc.proposal_textinput
+    proposal_textinput.on_change("name", proposal_textinput_callback)
+
+    def upload_button_callback(_attr, _old, new):
+        nonlocal cami_meta
+        with io.StringIO(base64.b64decode(new).decode()) as file:
+            cami_meta = pyzebra.parse_h5meta(file)
+        data_source.value = "cami file"
+        file_select_update()
+
+    upload_div = Div(text="or upload .cami file:", margin=(5, 5, 0, 5))
+    upload_button = FileInput(accept=".cami", width=200)
+    upload_button.on_change("value", upload_button_callback)
+
+    file_select = MultiSelect(title="Available .hdf files:", width=210, height=320)
+
+    def _init_datatable():
+        file_list = []
+        for scan in dataset:
+            file_list.append(os.path.basename(scan["original_filename"]))
+
+        scan_table_source.data.update(
+            file=file_list,
+            param=[None] * len(dataset),
+            frame=[None] * len(dataset),
+            x_pos=[None] * len(dataset),
+            y_pos=[None] * len(dataset),
+        )
+        scan_table_source.selected.indices = []
+        scan_table_source.selected.indices = [0]
+
+        param_select.value = "user defined"
+
+    def _update_table():
+        frame = []
+        x_pos = []
+        y_pos = []
+        for scan in dataset:
+            if "fit" in scan:
+                framei = scan["fit"]["frame"]
+                x_posi = scan["fit"]["x_pos"]
+                y_posi = scan["fit"]["y_pos"]
+            else:
+                framei = x_posi = y_posi = None
+
+            frame.append(framei)
+            x_pos.append(x_posi)
+            y_pos.append(y_posi)
+
+        scan_table_source.data.update(frame=frame, x_pos=x_pos, y_pos=y_pos)
+
+    def _file_open():
+        new_data = []
+        for f_name in file_select.value:
+            try:
+                new_data.append(pyzebra.read_detector_data(f_name))
+            except KeyError:
+                print("Could not read data from the file.")
+                return
+
+        dataset.extend(new_data)
+
+        _init_datatable()
+
+    def file_open_button_callback():
+        nonlocal dataset
+        dataset = []
+        _file_open()
+
+    file_open_button = Button(label="Open New", width=100)
+    file_open_button.on_click(file_open_button_callback)
+
+    def file_append_button_callback():
+        _file_open()
+
+    file_append_button = Button(label="Append", width=100)
+    file_append_button.on_click(file_append_button_callback)
+
+    # Scan select
+    def scan_table_select_callback(_attr, old, new):
+        if not new:
+            # skip empty selections
+            return
+
+        # Avoid selection of multiple indicies (via Shift+Click or Ctrl+Click)
+        if len(new) > 1:
+            # drop selection to the previous one
+            scan_table_source.selected.indices = old
+            return
+
+        if len(old) > 1:
+            # skip unnecessary update caused by selection drop
+            return
+
+        scan = dataset[new[0]]
+
+        zebra_mode = scan["zebra_mode"]
+        if zebra_mode == "nb":
+            metadata_table_source.data.update(geom=["normal beam"])
+        else:  # zebra_mode == "bi"
+            metadata_table_source.data.update(geom=["bisecting"])
+
+        if "mf" in scan:
+            metadata_table_source.data.update(mf=[scan["mf"][0]])
+        else:
+            metadata_table_source.data.update(mf=[None])
+
+        if "temp" in scan:
+            metadata_table_source.data.update(temp=[scan["temp"][0]])
+        else:
+            metadata_table_source.data.update(temp=[None])
+
+        update_overview_plot()
+
+    def scan_table_source_callback(_attr, _old, _new):
+        pass
+
+    scan_table_source = ColumnDataSource(dict(file=[], param=[], frame=[], x_pos=[], y_pos=[]))
+    scan_table_source.selected.on_change("indices", scan_table_select_callback)
+    scan_table_source.on_change("data", scan_table_source_callback)
+
+    scan_table = DataTable(
+        source=scan_table_source,
+        columns=[
+            TableColumn(field="file", title="file", editor=CellEditor(), width=150),
+            TableColumn(
+                field="param",
+                title="param",
+                formatter=num_formatter,
+                editor=NumberEditor(),
+                width=50,
+            ),
+            TableColumn(
+                field="frame", title="Frame", formatter=num_formatter, editor=CellEditor(), width=70
+            ),
+            TableColumn(
+                field="x_pos", title="X", formatter=num_formatter, editor=CellEditor(), width=70
+            ),
+            TableColumn(
+                field="y_pos", title="Y", formatter=num_formatter, editor=CellEditor(), width=70
+            ),
+        ],
+        width=470,  # +60 because of the index column
+        height=420,
+        editable=True,
+        autosize_mode="none",
+    )
+
+    def _get_selected_scan():
+        return dataset[scan_table_source.selected.indices[0]]
+
+    def param_select_callback(_attr, _old, new):
+        if new == "user defined":
+            param = [None] * len(dataset)
+        else:
+            # TODO: which value to take?
+            param = [scan[new][0] for scan in dataset]
+
+        scan_table_source.data["param"] = param
+        _update_param_plot()
+
+    param_select = Select(
+        title="Parameter:",
+        options=["user defined", "temp", "mf", "h", "k", "l"],
+        value="user defined",
+        width=145,
+    )
+    param_select.on_change("value", param_select_callback)
+
+    def update_overview_plot():
+        scan = _get_selected_scan()
+        counts = scan["counts"]
+        n_im, n_y, n_x = counts.shape
+        overview_x = np.mean(counts, axis=1)
+        overview_y = np.mean(counts, axis=2)
+
+        # normalize for simpler colormapping
+        overview_max_val = max(np.max(overview_x), np.max(overview_y))
+        overview_x = 1000 * overview_x / overview_max_val
+        overview_y = 1000 * overview_y / overview_max_val
+
+        overview_plot_x_image_source.data.update(image=[overview_x], dw=[n_x], dh=[n_im])
+        overview_plot_y_image_source.data.update(image=[overview_y], dw=[n_y], dh=[n_im])
+
+        if proj_auto_checkbox.active:
+            im_min = min(np.min(overview_x), np.min(overview_y))
+            im_max = max(np.max(overview_x), np.max(overview_y))
+
+            proj_display_min_spinner.value = im_min
+            proj_display_max_spinner.value = im_max
+
+            overview_plot_x_image_glyph.color_mapper.low = im_min
+            overview_plot_y_image_glyph.color_mapper.low = im_min
+            overview_plot_x_image_glyph.color_mapper.high = im_max
+            overview_plot_y_image_glyph.color_mapper.high = im_max
+
+        frame_range.start = 0
+        frame_range.end = n_im
+        frame_range.reset_start = 0
+        frame_range.reset_end = n_im
+        frame_range.bounds = (0, n_im)
+
+        scan_motor = scan["scan_motor"]
+        overview_plot_y.axis[1].axis_label = f"Scanning motor, {scan_motor}"
+
+        var = scan[scan_motor]
+        var_start = var[0]
+        var_end = var[-1] + (var[-1] - var[0]) / (n_im - 1)
+
+        scanning_motor_range.start = var_start
+        scanning_motor_range.end = var_end
+        scanning_motor_range.reset_start = var_start
+        scanning_motor_range.reset_end = var_end
+        # handle both, ascending and descending sequences
+        scanning_motor_range.bounds = (min(var_start, var_end), max(var_start, var_end))
+
+    # shared frame ranges
+    frame_range = Range1d(0, 1, bounds=(0, 1))
+    scanning_motor_range = Range1d(0, 1, bounds=(0, 1))
+
+    det_x_range = Range1d(0, IMAGE_W, bounds=(0, IMAGE_W))
+    overview_plot_x = Plot(
+        title=Title(text="Projections on X-axis"),
+        x_range=det_x_range,
+        y_range=frame_range,
+        extra_y_ranges={"scanning_motor": scanning_motor_range},
+        plot_height=400,
+        plot_width=IMAGE_PLOT_W - 3,
+    )
+
+    # ---- tools
+    wheelzoomtool = WheelZoomTool(maintain_focus=False)
+    overview_plot_x.toolbar.logo = None
+    overview_plot_x.add_tools(
+        PanTool(), BoxZoomTool(), wheelzoomtool, ResetTool(),
+    )
+    overview_plot_x.toolbar.active_scroll = wheelzoomtool
+
+    # ---- axes
+    overview_plot_x.add_layout(LinearAxis(axis_label="Coordinate X, pix"), place="below")
+    overview_plot_x.add_layout(
+        LinearAxis(axis_label="Frame", major_label_orientation="vertical"), place="left"
+    )
+
+    # ---- grid lines
+    overview_plot_x.add_layout(Grid(dimension=0, ticker=BasicTicker()))
+    overview_plot_x.add_layout(Grid(dimension=1, ticker=BasicTicker()))
+
+    # ---- rgba image glyph
+    overview_plot_x_image_source = ColumnDataSource(
+        dict(image=[np.zeros((1, 1), dtype="float32")], x=[0], y=[0], dw=[IMAGE_W], dh=[1])
+    )
+
+    overview_plot_x_image_glyph = Image(image="image", x="x", y="y", dw="dw", dh="dh")
+    overview_plot_x.add_glyph(
+        overview_plot_x_image_source, overview_plot_x_image_glyph, name="image_glyph"
+    )
+
+    det_y_range = Range1d(0, IMAGE_H, bounds=(0, IMAGE_H))
+    overview_plot_y = Plot(
+        title=Title(text="Projections on Y-axis"),
+        x_range=det_y_range,
+        y_range=frame_range,
+        extra_y_ranges={"scanning_motor": scanning_motor_range},
+        plot_height=400,
+        plot_width=IMAGE_PLOT_H + 22,
+    )
+
+    # ---- tools
+    wheelzoomtool = WheelZoomTool(maintain_focus=False)
+    overview_plot_y.toolbar.logo = None
+    overview_plot_y.add_tools(
+        PanTool(), BoxZoomTool(), wheelzoomtool, ResetTool(),
+    )
+    overview_plot_y.toolbar.active_scroll = wheelzoomtool
+
+    # ---- axes
+    overview_plot_y.add_layout(LinearAxis(axis_label="Coordinate Y, pix"), place="below")
+    overview_plot_y.add_layout(
+        LinearAxis(
+            y_range_name="scanning_motor",
+            axis_label="Scanning motor",
+            major_label_orientation="vertical",
+        ),
+        place="right",
+    )
+
+    # ---- grid lines
+    overview_plot_y.add_layout(Grid(dimension=0, ticker=BasicTicker()))
+    overview_plot_y.add_layout(Grid(dimension=1, ticker=BasicTicker()))
+
+    # ---- rgba image glyph
+    overview_plot_y_image_source = ColumnDataSource(
+        dict(image=[np.zeros((1, 1), dtype="float32")], x=[0], y=[0], dw=[IMAGE_H], dh=[1])
+    )
+
+    overview_plot_y_image_glyph = Image(image="image", x="x", y="y", dw="dw", dh="dh")
+    overview_plot_y.add_glyph(
+        overview_plot_y_image_source, overview_plot_y_image_glyph, name="image_glyph"
+    )
+
+    cmap_dict = {
+        "gray": Greys256,
+        "gray_reversed": Greys256[::-1],
+        "plasma": Plasma256,
+        "cividis": Cividis256,
+    }
+
+    def colormap_callback(_attr, _old, new):
+        overview_plot_x_image_glyph.color_mapper = LinearColorMapper(palette=cmap_dict[new])
+        overview_plot_y_image_glyph.color_mapper = LinearColorMapper(palette=cmap_dict[new])
+
+    colormap = Select(title="Colormap:", options=list(cmap_dict.keys()), width=210)
+    colormap.on_change("value", colormap_callback)
+    colormap.value = "plasma"
+
+    PROJ_STEP = 1
+
+    def proj_auto_checkbox_callback(state):
+        if state:
+            proj_display_min_spinner.disabled = True
+            proj_display_max_spinner.disabled = True
+        else:
+            proj_display_min_spinner.disabled = False
+            proj_display_max_spinner.disabled = False
+
+        update_overview_plot()
+
+    proj_auto_checkbox = CheckboxGroup(
+        labels=["Projections Intensity Range"], active=[0], width=145, margin=[10, 5, 0, 5]
+    )
+    proj_auto_checkbox.on_click(proj_auto_checkbox_callback)
+
+    def proj_display_max_spinner_callback(_attr, _old_value, new_value):
+        proj_display_min_spinner.high = new_value - PROJ_STEP
+        overview_plot_x_image_glyph.color_mapper.high = new_value
+        overview_plot_y_image_glyph.color_mapper.high = new_value
+
+    proj_display_max_spinner = Spinner(
+        low=0 + PROJ_STEP,
+        value=1,
+        step=PROJ_STEP,
+        disabled=bool(proj_auto_checkbox.active),
+        width=100,
+        height=31,
+    )
+    proj_display_max_spinner.on_change("value", proj_display_max_spinner_callback)
+
+    def proj_display_min_spinner_callback(_attr, _old_value, new_value):
+        proj_display_max_spinner.low = new_value + PROJ_STEP
+        overview_plot_x_image_glyph.color_mapper.low = new_value
+        overview_plot_y_image_glyph.color_mapper.low = new_value
+
+    proj_display_min_spinner = Spinner(
+        low=0,
+        high=1 - PROJ_STEP,
+        value=0,
+        step=PROJ_STEP,
+        disabled=bool(proj_auto_checkbox.active),
+        width=100,
+        height=31,
+    )
+    proj_display_min_spinner.on_change("value", proj_display_min_spinner_callback)
+
+    metadata_table_source = ColumnDataSource(dict(geom=[""], temp=[None], mf=[None]))
+    metadata_table = DataTable(
+        source=metadata_table_source,
+        columns=[
+            TableColumn(field="geom", title="Geometry", width=100),
+            TableColumn(field="temp", title="Temperature", formatter=num_formatter, width=100),
+            TableColumn(field="mf", title="Magnetic Field", formatter=num_formatter, width=100),
+        ],
+        width=300,
+        height=50,
+        autosize_mode="none",
+        index_position=None,
+    )
+
+    def _update_param_plot():
+        x = []
+        y = []
+        fit_param = fit_param_select.value
+        for s, p in zip(dataset, scan_table_source.data["param"]):
+            if "fit" in s and fit_param:
+                x.append(p)
+                y.append(s["fit"][fit_param])
+        param_plot_scatter_source.data.update(x=x, y=y)
+
+    # Parameter plot
+    param_plot = Plot(x_range=DataRange1d(), y_range=DataRange1d(), plot_height=400, plot_width=700)
+
+    param_plot.add_layout(LinearAxis(axis_label="Fit parameter"), place="left")
+    param_plot.add_layout(LinearAxis(axis_label="Parameter"), place="below")
+
+    param_plot.add_layout(Grid(dimension=0, ticker=BasicTicker()))
+    param_plot.add_layout(Grid(dimension=1, ticker=BasicTicker()))
+
+    param_plot_scatter_source = ColumnDataSource(dict(x=[], y=[]))
+    param_plot.add_glyph(param_plot_scatter_source, Scatter(x="x", y="y"))
+
+    param_plot.add_tools(PanTool(), WheelZoomTool(), ResetTool())
+    param_plot.toolbar.logo = None
+
+    def fit_param_select_callback(_attr, _old, _new):
+        _update_param_plot()
+
+    fit_param_select = Select(title="Fit parameter", options=[], width=145)
+    fit_param_select.on_change("value", fit_param_select_callback)
+
+    def proc_all_button_callback():
+        for scan in dataset:
+            pyzebra.fit_event(
+                scan,
+                int(np.floor(frame_range.start)),
+                int(np.ceil(frame_range.end)),
+                int(np.floor(det_y_range.start)),
+                int(np.ceil(det_y_range.end)),
+                int(np.floor(det_x_range.start)),
+                int(np.ceil(det_x_range.end)),
+            )
+
+        _update_table()
+
+        for scan in dataset:
+            if "fit" in scan:
+                options = list(scan["fit"].keys())
+                fit_param_select.options = options
+                fit_param_select.value = options[0]
+                break
+
+        _update_param_plot()
+
+    proc_all_button = Button(label="Process All", button_type="primary", width=145)
+    proc_all_button.on_click(proc_all_button_callback)
+
+    def proc_button_callback():
+        scan = _get_selected_scan()
+        pyzebra.fit_event(
+            scan,
+            int(np.floor(frame_range.start)),
+            int(np.ceil(frame_range.end)),
+            int(np.floor(det_y_range.start)),
+            int(np.ceil(det_y_range.end)),
+            int(np.floor(det_x_range.start)),
+            int(np.ceil(det_x_range.end)),
+        )
+
+        _update_table()
+
+        for scan in dataset:
+            if "fit" in scan:
+                options = list(scan["fit"].keys())
+                fit_param_select.options = options
+                fit_param_select.value = options[0]
+                break
+
+        _update_param_plot()
+
+    proc_button = Button(label="Process Current", width=145)
+    proc_button.on_click(proc_button_callback)
+
+    layout_controls = row(
+        colormap,
+        column(proj_auto_checkbox, row(proj_display_min_spinner, proj_display_max_spinner)),
+        proc_button,
+        proc_all_button,
+    )
+
+    layout_overview = column(
+        gridplot(
+            [[overview_plot_x, overview_plot_y]],
+            toolbar_options=dict(logo=None),
+            merge_tools=True,
+            toolbar_location="left",
+        ),
+        layout_controls,
+    )
+
+    # Plot tabs
+    plots = Tabs(
+        tabs=[
+            Panel(child=layout_overview, title="single scan"),
+            Panel(child=column(param_plot, row(fit_param_select)), title="parameter plot"),
+        ]
+    )
+
+    # Final layout
+    import_layout = column(
+        data_source,
+        upload_div,
+        upload_button,
+        file_select,
+        row(file_open_button, file_append_button),
+    )
+
+    scan_layout = column(scan_table, row(param_select, metadata_table))
+
+    tab_layout = column(row(import_layout, scan_layout, plots))
+
+    return Panel(child=tab_layout, title="hdf param study")

pyzebra/app/panel_param_study.py

@@ -0,0 +1,842 @@
+import base64
+import io
+import itertools
+import os
+import tempfile
+import types
+
+import numpy as np
+from bokeh.io import curdoc
+from bokeh.layouts import column, row
+from bokeh.models import (
+    BasicTicker,
+    Button,
+    CellEditor,
+    CheckboxEditor,
+    CheckboxGroup,
+    ColumnDataSource,
+    CustomJS,
+    DataRange1d,
+    DataTable,
+    Div,
+    Dropdown,
+    FileInput,
+    Grid,
+    HoverTool,
+    Image,
+    Legend,
+    Line,
+    LinearAxis,
+    LinearColorMapper,
+    MultiLine,
+    MultiSelect,
+    NumberEditor,
+    Panel,
+    PanTool,
+    Plot,
+    RadioGroup,
+    Range1d,
+    ResetTool,
+    Scatter,
+    Select,
+    Spacer,
+    Span,
+    Spinner,
+    TableColumn,
+    Tabs,
+    TextAreaInput,
+    WheelZoomTool,
+    Whisker,
+)
+from bokeh.palettes import Category10, Plasma256
+from scipy import interpolate
+
+import pyzebra
+from pyzebra.ccl_process import AREA_METHODS
+
+javaScript = """
+let j = 0;
+for (let i = 0; i < js_data.data['fname'].length; i++) {
+    if (js_data.data['content'][i] === "") continue;
+
+    setTimeout(function() {
+        const blob = new Blob([js_data.data['content'][i]], {type: 'text/plain'})
+        const link = document.createElement('a');
+        document.body.appendChild(link);
+        const url = window.URL.createObjectURL(blob);
+        link.href = url;
+        link.download = js_data.data['fname'][i] + js_data.data['ext'][i];
+        link.click();
+        window.URL.revokeObjectURL(url);
+        document.body.removeChild(link);
+    }, 100 * j)
+
+    j++;
+}
+"""
+
+
+def color_palette(n_colors):
+    palette = itertools.cycle(Category10[10])
+    return list(itertools.islice(palette, n_colors))
+
+
+def create():
+    doc = curdoc()
+    dataset = []
+    fit_params = {}
+    js_data = ColumnDataSource(data=dict(content=[""], fname=[""], ext=[""]))
+
+    def file_select_update_for_proposal():
+        proposal_path = proposal_textinput.name
+        if proposal_path:
+            file_list = []
+            for file in os.listdir(proposal_path):
+                if file.endswith((".ccl", ".dat")):
+                    file_list.append((os.path.join(proposal_path, file), file))
+            file_select.options = file_list
+            file_open_button.disabled = False
+            file_append_button.disabled = False
+        else:
+            file_select.options = []
+            file_open_button.disabled = True
+            file_append_button.disabled = True
+
+    doc.add_periodic_callback(file_select_update_for_proposal, 5000)
+
+    def proposal_textinput_callback(_attr, _old, _new):
+        file_select_update_for_proposal()
+
+    proposal_textinput = doc.proposal_textinput
+    proposal_textinput.on_change("name", proposal_textinput_callback)
+
+    def _init_datatable():
+        scan_list = [s["idx"] for s in dataset]
+        export = [s["export"] for s in dataset]
+        if param_select.value == "user defined":
+            param = [None] * len(dataset)
+        else:
+            param = [scan[param_select.value] for scan in dataset]
+
+        file_list = []
+        for scan in dataset:
+            file_list.append(os.path.basename(scan["original_filename"]))
+
+        scan_table_source.data.update(
+            file=file_list, scan=scan_list, param=param, fit=[0] * len(scan_list), export=export,
+        )
+        scan_table_source.selected.indices = []
+        scan_table_source.selected.indices = [0]
+
+        scan_motor_select.options = dataset[0]["scan_motors"]
+        scan_motor_select.value = dataset[0]["scan_motor"]
+
+        merge_options = [(str(i), f"{i} ({idx})") for i, idx in enumerate(scan_list)]
+        merge_from_select.options = merge_options
+        merge_from_select.value = merge_options[0][0]
+
+    file_select = MultiSelect(title="Available .ccl/.dat files:", width=210, height=250)
+
+    def file_open_button_callback():
+        nonlocal dataset
+        new_data = []
+        for f_path in file_select.value:
+            with open(f_path) as file:
+                f_name = os.path.basename(f_path)
+                base, ext = os.path.splitext(f_name)
+                try:
+                    file_data = pyzebra.parse_1D(file, ext)
+                except:
+                    print(f"Error loading {f_name}")
+                    continue
+
+            pyzebra.normalize_dataset(file_data, monitor_spinner.value)
+
+            if not new_data:  # first file
+                new_data = file_data
+                pyzebra.merge_duplicates(new_data)
+                js_data.data.update(fname=[base])
+            else:
+                pyzebra.merge_datasets(new_data, file_data)
+
+        if new_data:
+            dataset = new_data
+            _init_datatable()
+            append_upload_button.disabled = False
+
+    file_open_button = Button(label="Open New", width=100, disabled=True)
+    file_open_button.on_click(file_open_button_callback)
+
+    def file_append_button_callback():
+        file_data = []
+        for f_path in file_select.value:
+            with open(f_path) as file:
+                f_name = os.path.basename(f_path)
+                _, ext = os.path.splitext(f_name)
+                try:
+                    file_data = pyzebra.parse_1D(file, ext)
+                except:
+                    print(f"Error loading {f_name}")
+                    continue
+
+            pyzebra.normalize_dataset(file_data, monitor_spinner.value)
+            pyzebra.merge_datasets(dataset, file_data)
+
+        if file_data:
+            _init_datatable()
+
+    file_append_button = Button(label="Append", width=100, disabled=True)
+    file_append_button.on_click(file_append_button_callback)
+
+    def upload_button_callback(_attr, _old, _new):
+        nonlocal dataset
+        new_data = []
+        for f_str, f_name in zip(upload_button.value, upload_button.filename):
+            with io.StringIO(base64.b64decode(f_str).decode()) as file:
+                base, ext = os.path.splitext(f_name)
+                try:
+                    file_data = pyzebra.parse_1D(file, ext)
+                except:
+                    print(f"Error loading {f_name}")
+                    continue
+
+            pyzebra.normalize_dataset(file_data, monitor_spinner.value)
+
+            if not new_data:  # first file
+                new_data = file_data
+                pyzebra.merge_duplicates(new_data)
+                js_data.data.update(fname=[base])
+            else:
+                pyzebra.merge_datasets(new_data, file_data)
+
+        if new_data:
+            dataset = new_data
+            _init_datatable()
+            append_upload_button.disabled = False
+
+    upload_div = Div(text="or upload new .ccl/.dat files:", margin=(5, 5, 0, 5))
+    upload_button = FileInput(accept=".ccl,.dat", multiple=True, width=200)
+    # for on_change("value", ...) or on_change("filename", ...),
+    # see https://github.com/bokeh/bokeh/issues/11461
+    upload_button.on_change("filename", upload_button_callback)
+
+    def append_upload_button_callback(_attr, _old, _new):
+        file_data = []
+        for f_str, f_name in zip(append_upload_button.value, append_upload_button.filename):
+            with io.StringIO(base64.b64decode(f_str).decode()) as file:
+                _, ext = os.path.splitext(f_name)
+                try:
+                    file_data = pyzebra.parse_1D(file, ext)
+                except:
+                    print(f"Error loading {f_name}")
+                    continue
+
+            pyzebra.normalize_dataset(file_data, monitor_spinner.value)
+            pyzebra.merge_datasets(dataset, file_data)
+
+        if file_data:
+            _init_datatable()
+
+    append_upload_div = Div(text="append extra files:", margin=(5, 5, 0, 5))
+    append_upload_button = FileInput(accept=".ccl,.dat", multiple=True, width=200, disabled=True)
+    # for on_change("value", ...) or on_change("filename", ...),
+    # see https://github.com/bokeh/bokeh/issues/11461
+    append_upload_button.on_change("filename", append_upload_button_callback)
+
+    def monitor_spinner_callback(_attr, _old, new):
+        if dataset:
+            pyzebra.normalize_dataset(dataset, new)
+            _update_single_scan_plot()
+            _update_overview()
+
+    monitor_spinner = Spinner(title="Monitor:", mode="int", value=100_000, low=1, width=145)
+    monitor_spinner.on_change("value", monitor_spinner_callback)
+
+    def scan_motor_select_callback(_attr, _old, new):
+        if dataset:
+            for scan in dataset:
+                scan["scan_motor"] = new
+            _update_single_scan_plot()
+            _update_overview()
+
+    scan_motor_select = Select(title="Scan motor:", options=[], width=145)
+    scan_motor_select.on_change("value", scan_motor_select_callback)
+
+    def _update_table():
+        fit_ok = [(1 if "fit" in scan else 0) for scan in dataset]
+        export = [scan["export"] for scan in dataset]
+        if param_select.value == "user defined":
+            param = [None] * len(dataset)
+        else:
+            param = [scan[param_select.value] for scan in dataset]
+
+        scan_table_source.data.update(fit=fit_ok, export=export, param=param)
+
+    def _update_single_scan_plot():
+        scan = _get_selected_scan()
+        scan_motor = scan["scan_motor"]
+
+        y = scan["counts"]
+        y_err = scan["counts_err"]
+        x = scan[scan_motor]
+
+        plot.axis[0].axis_label = scan_motor
+        plot_scatter_source.data.update(x=x, y=y, y_upper=y + y_err, y_lower=y - y_err)
+
+        fit = scan.get("fit")
+        if fit is not None:
+            x_fit = np.linspace(x[0], x[-1], 100)
+            plot_fit_source.data.update(x=x_fit, y=fit.eval(x=x_fit))
+
+            x_bkg = []
+            y_bkg = []
+            xs_peak = []
+            ys_peak = []
+            comps = fit.eval_components(x=x_fit)
+            for i, model in enumerate(fit_params):
+                if "linear" in model:
+                    x_bkg = x_fit
+                    y_bkg = comps[f"f{i}_"]
+
+                elif any(val in model for val in ("gaussian", "voigt", "pvoigt")):
+                    xs_peak.append(x_fit)
+                    ys_peak.append(comps[f"f{i}_"])
+
+            plot_bkg_source.data.update(x=x_bkg, y=y_bkg)
+            plot_peak_source.data.update(xs=xs_peak, ys=ys_peak)
+
+            fit_output_textinput.value = fit.fit_report()
+
+        else:
+            plot_fit_source.data.update(x=[], y=[])
+            plot_bkg_source.data.update(x=[], y=[])
+            plot_peak_source.data.update(xs=[], ys=[])
+            fit_output_textinput.value = ""
+
+    def _update_overview():
+        xs = []
+        ys = []
+        param = []
+        x = []
+        y = []
+        par = []
+        for s, p in enumerate(scan_table_source.data["param"]):
+            if p is not None:
+                scan = dataset[s]
+                scan_motor = scan["scan_motor"]
+                xs.append(scan[scan_motor])
+                x.extend(scan[scan_motor])
+                ys.append(scan["counts"])
+                y.extend([float(p)] * len(scan[scan_motor]))
+                param.append(float(p))
+                par.extend(scan["counts"])
+
+        if dataset:
+            scan_motor = dataset[0]["scan_motor"]
+            ov_plot.axis[0].axis_label = scan_motor
+            ov_param_plot.axis[0].axis_label = scan_motor
+
+        ov_plot_mline_source.data.update(xs=xs, ys=ys, param=param, color=color_palette(len(xs)))
+
+        ov_param_plot_scatter_source.data.update(x=x, y=y)
+
+        if y:
+            x1, x2 = min(x), max(x)
+            y1, y2 = min(y), max(y)
+            grid_x, grid_y = np.meshgrid(
+                np.linspace(x1, x2, ov_param_plot.inner_width),
+                np.linspace(y1, y2, ov_param_plot.inner_height),
+            )
+            image = interpolate.griddata((x, y), par, (grid_x, grid_y))
+            ov_param_plot_image_source.data.update(
+                image=[image], x=[x1], y=[y1], dw=[x2 - x1], dh=[y2 - y1]
+            )
+
+            x_range = ov_param_plot.x_range
+            x_range.start, x_range.end = x1, x2
+            x_range.reset_start, x_range.reset_end = x1, x2
+            x_range.bounds = (x1, x2)
+
+            y_range = ov_param_plot.y_range
+            y_range.start, y_range.end = y1, y2
+            y_range.reset_start, y_range.reset_end = y1, y2
+            y_range.bounds = (y1, y2)
+
+        else:
+            ov_param_plot_image_source.data.update(image=[], x=[], y=[], dw=[], dh=[])
+
+    def _update_param_plot():
+        x = []
+        y = []
+        y_lower = []
+        y_upper = []
+        fit_param = fit_param_select.value
+        for s, p in zip(dataset, scan_table_source.data["param"]):
+            if "fit" in s and fit_param:
+                x.append(p)
+                param_fit_val = s["fit"].params[fit_param].value
+                param_fit_std = s["fit"].params[fit_param].stderr
+                if param_fit_std is None:
+                    param_fit_std = 0
+                y.append(param_fit_val)
+                y_lower.append(param_fit_val - param_fit_std)
+                y_upper.append(param_fit_val + param_fit_std)
+
+        param_plot_scatter_source.data.update(x=x, y=y, y_lower=y_lower, y_upper=y_upper)
+
+    # Main plot
+    plot = Plot(
+        x_range=DataRange1d(),
+        y_range=DataRange1d(only_visible=True),
+        plot_height=450,
+        plot_width=700,
+    )
+
+    plot.add_layout(LinearAxis(axis_label="Counts"), place="left")
+    plot.add_layout(LinearAxis(axis_label="Scan motor"), place="below")
+
+    plot.add_layout(Grid(dimension=0, ticker=BasicTicker()))
+    plot.add_layout(Grid(dimension=1, ticker=BasicTicker()))
+
+    plot_scatter_source = ColumnDataSource(dict(x=[0], y=[0], y_upper=[0], y_lower=[0]))
+    plot_scatter = plot.add_glyph(
+        plot_scatter_source, Scatter(x="x", y="y", line_color="steelblue", fill_color="steelblue")
+    )
+    plot.add_layout(Whisker(source=plot_scatter_source, base="x", upper="y_upper", lower="y_lower"))
+
+    plot_fit_source = ColumnDataSource(dict(x=[0], y=[0]))
+    plot_fit = plot.add_glyph(plot_fit_source, Line(x="x", y="y"))
+
+    plot_bkg_source = ColumnDataSource(dict(x=[0], y=[0]))
+    plot_bkg = plot.add_glyph(
+        plot_bkg_source, Line(x="x", y="y", line_color="green", line_dash="dashed")
+    )
+
+    plot_peak_source = ColumnDataSource(dict(xs=[[0]], ys=[[0]]))
+    plot_peak = plot.add_glyph(
+        plot_peak_source, MultiLine(xs="xs", ys="ys", line_color="red", line_dash="dashed")
+    )
+
+    fit_from_span = Span(location=None, dimension="height", line_dash="dashed")
+    plot.add_layout(fit_from_span)
+
+    fit_to_span = Span(location=None, dimension="height", line_dash="dashed")
+    plot.add_layout(fit_to_span)
+
+    plot.add_layout(
+        Legend(
+            items=[
+                ("data", [plot_scatter]),
+                ("best fit", [plot_fit]),
+                ("peak", [plot_peak]),
+                ("linear", [plot_bkg]),
+            ],
+            location="top_left",
+            click_policy="hide",
+        )
+    )
+
+    plot.add_tools(PanTool(), WheelZoomTool(), ResetTool())
+    plot.toolbar.logo = None
+
+    # Overview multilines plot
+    ov_plot = Plot(x_range=DataRange1d(), y_range=DataRange1d(), plot_height=450, plot_width=700)
+
+    ov_plot.add_layout(LinearAxis(axis_label="Counts"), place="left")
+    ov_plot.add_layout(LinearAxis(axis_label="Scan motor"), place="below")
+
+    ov_plot.add_layout(Grid(dimension=0, ticker=BasicTicker()))
+    ov_plot.add_layout(Grid(dimension=1, ticker=BasicTicker()))
+
+    ov_plot_mline_source = ColumnDataSource(dict(xs=[], ys=[], param=[], color=[]))
+    ov_plot.add_glyph(ov_plot_mline_source, MultiLine(xs="xs", ys="ys", line_color="color"))
+
+    hover_tool = HoverTool(tooltips=[("param", "@param")])
+    ov_plot.add_tools(PanTool(), WheelZoomTool(), hover_tool, ResetTool())
+
+    ov_plot.add_tools(PanTool(), WheelZoomTool(), ResetTool())
+    ov_plot.toolbar.logo = None
+
+    # Overview perams plot
+    ov_param_plot = Plot(x_range=Range1d(), y_range=Range1d(), plot_height=450, plot_width=700)
+
+    ov_param_plot.add_layout(LinearAxis(axis_label="Param"), place="left")
+    ov_param_plot.add_layout(LinearAxis(axis_label="Scan motor"), place="below")
+
+    ov_param_plot.add_layout(Grid(dimension=0, ticker=BasicTicker()))
+    ov_param_plot.add_layout(Grid(dimension=1, ticker=BasicTicker()))
+
+    color_mapper = LinearColorMapper(palette=Plasma256)
+    ov_param_plot_image_source = ColumnDataSource(dict(image=[], x=[], y=[], dw=[], dh=[]))
+    ov_param_plot.add_glyph(
+        ov_param_plot_image_source,
+        Image(image="image", x="x", y="y", dw="dw", dh="dh", color_mapper=color_mapper),
+    )
+
+    ov_param_plot_scatter_source = ColumnDataSource(dict(x=[], y=[]))
+    ov_param_plot.add_glyph(
+        ov_param_plot_scatter_source, Scatter(x="x", y="y", marker="dot", size=15),
+    )
+
+    ov_param_plot.add_tools(PanTool(), WheelZoomTool(), ResetTool())
+    ov_param_plot.toolbar.logo = None
+
+    # Parameter plot
+    param_plot = Plot(x_range=DataRange1d(), y_range=DataRange1d(), plot_height=400, plot_width=700)
+
+    param_plot.add_layout(LinearAxis(axis_label="Fit parameter"), place="left")
+    param_plot.add_layout(LinearAxis(axis_label="Parameter"), place="below")
+
+    param_plot.add_layout(Grid(dimension=0, ticker=BasicTicker()))
+    param_plot.add_layout(Grid(dimension=1, ticker=BasicTicker()))
+
+    param_plot_scatter_source = ColumnDataSource(dict(x=[], y=[], y_upper=[], y_lower=[]))
+    param_plot.add_glyph(param_plot_scatter_source, Scatter(x="x", y="y"))
+    param_plot.add_layout(
+        Whisker(source=param_plot_scatter_source, base="x", upper="y_upper", lower="y_lower")
+    )
+
+    param_plot.add_tools(PanTool(), WheelZoomTool(), ResetTool())
+    param_plot.toolbar.logo = None
+
+    def fit_param_select_callback(_attr, _old, _new):
+        _update_param_plot()
+
+    fit_param_select = Select(title="Fit parameter", options=[], width=145)
+    fit_param_select.on_change("value", fit_param_select_callback)
+
+    # Plot tabs
+    plots = Tabs(
+        tabs=[
+            Panel(child=plot, title="single scan"),
+            Panel(child=ov_plot, title="overview"),
+            Panel(child=ov_param_plot, title="overview map"),
+            Panel(child=column(param_plot, row(fit_param_select)), title="parameter plot"),
+        ]
+    )
+
+    # Scan select
+    def scan_table_select_callback(_attr, old, new):
+        if not new:
+            # skip empty selections
+            return
+
+        # Avoid selection of multiple indicies (via Shift+Click or Ctrl+Click)
+        if len(new) > 1:
+            # drop selection to the previous one
+            scan_table_source.selected.indices = old
+            return
+
+        if len(old) > 1:
+            # skip unnecessary update caused by selection drop
+            return
+
+        _update_single_scan_plot()
+
+    def scan_table_source_callback(_attr, _old, new):
+        # unfortunately, we don't know if the change comes from data update or user input
+        # also `old` and `new` are the same for non-scalars
+        for scan, export in zip(dataset, new["export"]):
+            scan["export"] = export
+        _update_overview()
+        _update_param_plot()
+        _update_preview()
+
+    scan_table_source = ColumnDataSource(dict(file=[], scan=[], param=[], fit=[], export=[]))
+    scan_table_source.on_change("data", scan_table_source_callback)
+    scan_table_source.selected.on_change("indices", scan_table_select_callback)
+
+    scan_table = DataTable(
+        source=scan_table_source,
+        columns=[
+            TableColumn(field="file", title="file", editor=CellEditor(), width=150),
+            TableColumn(field="scan", title="scan", editor=CellEditor(), width=50),
+            TableColumn(field="param", title="param", editor=NumberEditor(), width=50),
+            TableColumn(field="fit", title="Fit", editor=CellEditor(), width=50),
+            TableColumn(field="export", title="Export", editor=CheckboxEditor(), width=50),
+        ],
+        width=410,  # +60 because of the index column
+        height=350,
+        editable=True,
+        autosize_mode="none",
+    )
+
+    merge_from_select = Select(title="scan:", width=145)
+
+    def merge_button_callback():
+        scan_into = _get_selected_scan()
+        scan_from = dataset[int(merge_from_select.value)]
+
+        if scan_into is scan_from:
+            print("WARNING: Selected scans for merging are identical")
+            return
+
+        pyzebra.merge_scans(scan_into, scan_from)
+        _update_table()
+        _update_single_scan_plot()
+        _update_overview()
+
+    merge_button = Button(label="Merge into current", width=145)
+    merge_button.on_click(merge_button_callback)
+
+    def restore_button_callback():
+        pyzebra.restore_scan(_get_selected_scan())
+        _update_table()
+        _update_single_scan_plot()
+        _update_overview()
+
+    restore_button = Button(label="Restore scan", width=145)
+    restore_button.on_click(restore_button_callback)
+
+    def _get_selected_scan():
+        return dataset[scan_table_source.selected.indices[0]]
+
+    def param_select_callback(_attr, _old, _new):
+        _update_table()
+
+    param_select = Select(
+        title="Parameter:",
+        options=["user defined", "temp", "mf", "h", "k", "l"],
+        value="user defined",
+        width=145,
+    )
+    param_select.on_change("value", param_select_callback)
+
+    def fit_from_spinner_callback(_attr, _old, new):
+        fit_from_span.location = new
+
+    fit_from_spinner = Spinner(title="Fit from:", width=145)
+    fit_from_spinner.on_change("value", fit_from_spinner_callback)
+
+    def fit_to_spinner_callback(_attr, _old, new):
+        fit_to_span.location = new
+
+    fit_to_spinner = Spinner(title="to:", width=145)
+    fit_to_spinner.on_change("value", fit_to_spinner_callback)
+
+    def fitparams_add_dropdown_callback(click):
+        # bokeh requires (str, str) for MultiSelect options
+        new_tag = f"{click.item}-{fitparams_select.tags[0]}"
+        fitparams_select.options.append((new_tag, click.item))
+        fit_params[new_tag] = fitparams_factory(click.item)
+        fitparams_select.tags[0] += 1
+
+    fitparams_add_dropdown = Dropdown(
+        label="Add fit function",
+        menu=[
+            ("Linear", "linear"),
+            ("Gaussian", "gaussian"),
+            ("Voigt", "voigt"),
+            ("Pseudo Voigt", "pvoigt"),
+            # ("Pseudo Voigt1", "pseudovoigt1"),
+        ],
+        width=145,
+    )
+    fitparams_add_dropdown.on_click(fitparams_add_dropdown_callback)
+
+    def fitparams_select_callback(_attr, old, new):
+        # Avoid selection of multiple indicies (via Shift+Click or Ctrl+Click)
+        if len(new) > 1:
+            # drop selection to the previous one
+            fitparams_select.value = old
+            return
+
+        if len(old) > 1:
+            # skip unnecessary update caused by selection drop
+            return
+
+        if new:
+            fitparams_table_source.data.update(fit_params[new[0]])
+        else:
+            fitparams_table_source.data.update(dict(param=[], value=[], vary=[], min=[], max=[]))
+
+    fitparams_select = MultiSelect(options=[], height=120, width=145)
+    fitparams_select.tags = [0]
+    fitparams_select.on_change("value", fitparams_select_callback)
+
+    def fitparams_remove_button_callback():
+        if fitparams_select.value:
+            sel_tag = fitparams_select.value[0]
+            del fit_params[sel_tag]
+            for elem in fitparams_select.options:
+                if elem[0] == sel_tag:
+                    fitparams_select.options.remove(elem)
+                    break
+
+            fitparams_select.value = []
+
+    fitparams_remove_button = Button(label="Remove fit function", width=145)
+    fitparams_remove_button.on_click(fitparams_remove_button_callback)
+
+    def fitparams_factory(function):
+        if function == "linear":
+            params = ["slope", "intercept"]
+        elif function == "gaussian":
+            params = ["amplitude", "center", "sigma"]
+        elif function == "voigt":
+            params = ["amplitude", "center", "sigma", "gamma"]
+        elif function == "pvoigt":
+            params = ["amplitude", "center", "sigma", "fraction"]
+        elif function == "pseudovoigt1":
+            params = ["amplitude", "center", "g_sigma", "l_sigma", "fraction"]
+        else:
+            raise ValueError("Unknown fit function")
+
+        n = len(params)
+        fitparams = dict(
+            param=params, value=[None] * n, vary=[True] * n, min=[None] * n, max=[None] * n,
+        )
+
+        if function == "linear":
+            fitparams["value"] = [0, 1]
+            fitparams["vary"] = [False, True]
+            fitparams["min"] = [None, 0]
+
+        elif function == "gaussian":
+            fitparams["min"] = [0, None, None]
+
+        return fitparams
+
+    fitparams_table_source = ColumnDataSource(dict(param=[], value=[], vary=[], min=[], max=[]))
+    fitparams_table = DataTable(
+        source=fitparams_table_source,
+        columns=[
+            TableColumn(field="param", title="Parameter", editor=CellEditor()),
+            TableColumn(field="value", title="Value", editor=NumberEditor()),
+            TableColumn(field="vary", title="Vary", editor=CheckboxEditor()),
+            TableColumn(field="min", title="Min", editor=NumberEditor()),
+            TableColumn(field="max", title="Max", editor=NumberEditor()),
+        ],
+        height=200,
+        width=350,
+        index_position=None,
+        editable=True,
+        auto_edit=True,
+    )
+
+    # start with `background` and `gauss` fit functions added
+    fitparams_add_dropdown_callback(types.SimpleNamespace(item="linear"))
+    fitparams_add_dropdown_callback(types.SimpleNamespace(item="gaussian"))
+    fitparams_select.value = ["gaussian-1"]  # add selection to gauss
+
+    fit_output_textinput = TextAreaInput(title="Fit results:", width=750, height=200)
+
+    def proc_all_button_callback():
+        for scan in dataset:
+            if scan["export"]:
+                pyzebra.fit_scan(
+                    scan, fit_params, fit_from=fit_from_spinner.value, fit_to=fit_to_spinner.value
+                )
+                pyzebra.get_area(
+                    scan,
+                    area_method=AREA_METHODS[area_method_radiobutton.active],
+                    lorentz=lorentz_checkbox.active,
+                )
+
+        _update_single_scan_plot()
+        _update_overview()
+        _update_table()
+
+        for scan in dataset:
+            if "fit" in scan:
+                options = list(scan["fit"].params.keys())
+                fit_param_select.options = options
+                fit_param_select.value = options[0]
+                break
+
+    proc_all_button = Button(label="Process All", button_type="primary", width=145)
+    proc_all_button.on_click(proc_all_button_callback)
+
+    def proc_button_callback():
+        scan = _get_selected_scan()
+        pyzebra.fit_scan(
+            scan, fit_params, fit_from=fit_from_spinner.value, fit_to=fit_to_spinner.value
+        )
+        pyzebra.get_area(
+            scan,
+            area_method=AREA_METHODS[area_method_radiobutton.active],
+            lorentz=lorentz_checkbox.active,
+        )
+
+        _update_single_scan_plot()
+        _update_overview()
+        _update_table()
+
+        for scan in dataset:
+            if "fit" in scan:
+                options = list(scan["fit"].params.keys())
+                fit_param_select.options = options
+                fit_param_select.value = options[0]
+                break
+
+    proc_button = Button(label="Process Current", width=145)
+    proc_button.on_click(proc_button_callback)
+
+    area_method_div = Div(text="Intensity:", margin=(5, 5, 0, 5))
+    area_method_radiobutton = RadioGroup(labels=["Function", "Area"], active=0, width=145)
+
+    lorentz_checkbox = CheckboxGroup(labels=["Lorentz Correction"], width=145, margin=(13, 5, 5, 5))
+
+    export_preview_textinput = TextAreaInput(title="Export file preview:", width=450, height=400)
+
+    def _update_preview():
+        with tempfile.TemporaryDirectory() as temp_dir:
+            temp_file = temp_dir + "/temp"
+            export_data = []
+            param_data = []
+            for scan, param in zip(dataset, scan_table_source.data["param"]):
+                if scan["export"] and param:
+                    export_data.append(scan)
+                    param_data.append(param)
+
+            pyzebra.export_param_study(export_data, param_data, temp_file)
+
+            exported_content = ""
+            file_content = []
+
+            fname = temp_file
+            if os.path.isfile(fname):
+                with open(fname) as f:
+                    content = f.read()
+                    exported_content += content
+            else:
+                content = ""
+            file_content.append(content)
+
+            js_data.data.update(content=file_content)
+            export_preview_textinput.value = exported_content
+
+    save_button = Button(label="Download File", button_type="success", width=220)
+    save_button.js_on_click(CustomJS(args={"js_data": js_data}, code=javaScript))
+
+    fitpeak_controls = row(
+        column(fitparams_add_dropdown, fitparams_select, fitparams_remove_button),
+        fitparams_table,
+        Spacer(width=20),
+        column(fit_from_spinner, lorentz_checkbox, area_method_div, area_method_radiobutton),
+        column(fit_to_spinner, proc_button, proc_all_button),
+    )
+
+    scan_layout = column(
+        scan_table,
+        row(monitor_spinner, scan_motor_select, param_select),
+        row(column(Spacer(height=19), row(restore_button, merge_button)), merge_from_select),
+    )
+
+    import_layout = column(
+        file_select,
+        row(file_open_button, file_append_button),
+        upload_div,
+        upload_button,
+        append_upload_div,
+        append_upload_button,
+    )
+
+    export_layout = column(export_preview_textinput, row(save_button))
+
+    tab_layout = column(
+        row(import_layout, scan_layout, plots, Spacer(width=30), export_layout),
+        row(fitpeak_controls, fit_output_textinput),
+    )
+
+    return Panel(child=tab_layout, title="param study")

pyzebra/app/panel_spind.py

@@ -0,0 +1,223 @@
+import os
+import subprocess
+import tempfile
+
+import numpy as np
+from bokeh.io import curdoc
+from bokeh.layouts import column, row
+from bokeh.models import (
+    Button,
+    ColumnDataSource,
+    DataTable,
+    Panel,
+    Spinner,
+    TableColumn,
+    TextAreaInput,
+    TextInput,
+)
+
+import pyzebra
+
+
+def create():
+    doc = curdoc()
+    events_data = doc.events_data
+
+    npeaks_spinner = Spinner(title="Number of peaks from hdf_view panel:", disabled=True)
+    lattice_const_textinput = TextInput(title="Lattice constants:")
+    max_res_spinner = Spinner(title="max-res:", value=2, step=0.01, width=145)
+    seed_pool_size_spinner = Spinner(title="seed-pool-size:", value=5, step=0.01, width=145)
+    seed_len_tol_spinner = Spinner(title="seed-len-tol:", value=0.02, step=0.01, width=145)
+    seed_angle_tol_spinner = Spinner(title="seed-angle-tol:", value=1, step=0.01, width=145)
+    eval_hkl_tol_spinner = Spinner(title="eval-hkl-tol:", value=0.15, step=0.01, width=145)
+
+    diff_vec = []
+    ub_matrices = []
+
+    def process_button_callback():
+        # drop table selection to clear result fields
+        results_table_source.selected.indices = []
+
+        nonlocal diff_vec
+        with tempfile.TemporaryDirectory() as temp_dir:
+            temp_peak_list_dir = os.path.join(temp_dir, "peak_list")
+            os.mkdir(temp_peak_list_dir)
+            temp_event_file = os.path.join(temp_peak_list_dir, "event-0.txt")
+            temp_hkl_file = os.path.join(temp_dir, "hkl.h5")
+
+            comp_proc = subprocess.run(
+                [
+                    "mpiexec",
+                    "-n",
+                    "2",
+                    "python",
+                    os.path.join(doc.spind_path, "gen_hkl_table.py"),
+                    lattice_const_textinput.value,
+                    "--max-res",
+                    str(max_res_spinner.value),
+                    "-o",
+                    temp_hkl_file,
+                ],
+                check=True,
+                stdout=subprocess.PIPE,
+                stderr=subprocess.STDOUT,
+                text=True,
+            )
+            print(" ".join(comp_proc.args))
+            print(comp_proc.stdout)
+
+            # prepare an event file
+            diff_vec = []
+            with open(temp_event_file, "w") as f:
+                npeaks = len(next(iter(doc.events_data.values())))
+                for ind in range(npeaks):
+                    wave = events_data["wave"][ind]
+                    ddist = events_data["ddist"][ind]
+                    x_pos = events_data["x_pos"][ind]
+                    y_pos = events_data["y_pos"][ind]
+                    intensity = events_data["intensity"][ind]
+                    snr_cnts = events_data["snr_cnts"][ind]
+                    gamma = events_data["gamma"][ind]
+                    omega = events_data["omega"][ind]
+                    chi = events_data["chi"][ind]
+                    phi = events_data["phi"][ind]
+                    nu = events_data["nu"][ind]
+
+                    ga, nu = pyzebra.det2pol(ddist, gamma, nu, x_pos, y_pos)
+                    diff_vector = pyzebra.z1frmd(wave, ga, omega, chi, phi, nu)
+                    d_spacing = float(pyzebra.dandth(wave, diff_vector)[0])
+                    diff_vector = diff_vector.flatten() * 1e10
+                    dv1, dv2, dv3 = diff_vector
+
+                    diff_vec.append(diff_vector)
+                    f.write(
+                        f"{x_pos} {y_pos} {intensity} {snr_cnts} {dv1} {dv2} {dv3} {d_spacing}\n"
+                    )
+
+            print(f"Content of {temp_event_file}:")
+            with open(temp_event_file) as f:
+                print(f.read())
+
+            comp_proc = subprocess.run(
+                [
+                    "mpiexec",
+                    "-n",
+                    "2",
+                    "python",
+                    os.path.join(doc.spind_path, "SPIND.py"),
+                    temp_peak_list_dir,
+                    temp_hkl_file,
+                    "-o",
+                    temp_dir,
+                    "--seed-pool-size",
+                    str(seed_pool_size_spinner.value),
+                    "--seed-len-tol",
+                    str(seed_len_tol_spinner.value),
+                    "--seed-angle-tol",
+                    str(seed_angle_tol_spinner.value),
+                    "--eval-hkl-tol",
+                    str(eval_hkl_tol_spinner.value),
+                ],
+                check=True,
+                stdout=subprocess.PIPE,
+                stderr=subprocess.STDOUT,
+                text=True,
+            )
+            print(" ".join(comp_proc.args))
+            print(comp_proc.stdout)
+
+            spind_out_file = os.path.join(temp_dir, "spind.txt")
+            spind_res = dict(
+                label=[], crystal_id=[], match_rate=[], matched_peaks=[], column_5=[], ub_matrix=[],
+            )
+            try:
+                with open(spind_out_file) as f_out:
+                    for line in f_out:
+                        c1, c2, c3, c4, c5, *c_rest = line.split()
+                        spind_res["label"].append(c1)
+                        spind_res["crystal_id"].append(c2)
+                        spind_res["match_rate"].append(c3)
+                        spind_res["matched_peaks"].append(c4)
+                        spind_res["column_5"].append(c5)
+
+                        # last digits are spind UB matrix
+                        vals = list(map(float, c_rest))
+                        ub_matrix_spind = np.transpose(np.array(vals).reshape(3, 3))
+                        ub_matrix = np.linalg.inv(ub_matrix_spind)
+                        ub_matrices.append(ub_matrix)
+                        spind_res["ub_matrix"].append(str(ub_matrix_spind * 1e-10))
+
+                print(f"Content of {spind_out_file}:")
+                with open(spind_out_file) as f:
+                    print(f.read())
+
+            except FileNotFoundError:
+                print("No results from spind")
+
+            results_table_source.data.update(spind_res)
+
+    process_button = Button(label="Process", button_type="primary")
+    process_button.on_click(process_button_callback)
+
+    if doc.spind_path is None:
+        process_button.disabled = True
+
+    ub_matrix_textareainput = TextAreaInput(title="UB matrix:", rows=7, width=400)
+    hkl_textareainput = TextAreaInput(title="hkl values:", rows=7, width=400)
+
+    def results_table_select_callback(_attr, old, new):
+        if new:
+            ind = new[0]
+            ub_matrix = ub_matrices[ind]
+            res = ""
+            for vec in diff_vec:
+                res += f"{ub_matrix @ vec}\n"
+            ub_matrix_textareainput.value = str(ub_matrix * 1e10)
+            hkl_textareainput.value = res
+        else:
+            ub_matrix_textareainput.value = ""
+            hkl_textareainput.value = ""
+
+    results_table_source = ColumnDataSource(
+        dict(label=[], crystal_id=[], match_rate=[], matched_peaks=[], column_5=[], ub_matrix=[])
+    )
+    results_table = DataTable(
+        source=results_table_source,
+        columns=[
+            TableColumn(field="label", title="Label", width=50),
+            TableColumn(field="crystal_id", title="Crystal ID", width=100),
+            TableColumn(field="match_rate", title="Match Rate", width=100),
+            TableColumn(field="matched_peaks", title="Matched Peaks", width=100),
+            TableColumn(field="column_5", title="", width=100),
+            TableColumn(field="ub_matrix", title="UB Matrix", width=700),
+        ],
+        height=300,
+        width=1200,
+        autosize_mode="none",
+        index_position=None,
+    )
+
+    results_table_source.selected.on_change("indices", results_table_select_callback)
+
+    tab_layout = row(
+        column(
+            npeaks_spinner,
+            lattice_const_textinput,
+            row(max_res_spinner, seed_pool_size_spinner),
+            row(seed_len_tol_spinner, seed_angle_tol_spinner),
+            row(eval_hkl_tol_spinner),
+            process_button,
+        ),
+        column(results_table, row(ub_matrix_textareainput, hkl_textareainput)),
+    )
+
+    async def update_npeaks_spinner():
+        npeaks = len(next(iter(doc.events_data.values())))
+        npeaks_spinner.value = npeaks
+        # TODO: check cell parameter for consistency?
+        if npeaks:
+            lattice_const_textinput.value = ",".join(map(str, doc.events_data["cell"][0]))
+
+    doc.add_periodic_callback(update_npeaks_spinner, 1000)
+
+    return Panel(child=tab_layout, title="spind")

pyzebra/ccl_findpeaks.py

@@ -1,81 +0,0 @@
-import numpy as np
-import scipy as sc
-from scipy.interpolate import interp1d
-from scipy.signal import savgol_filter
-
-
-def ccl_findpeaks(
-    scan,
-    int_threshold=0.8,
-    prominence=50,
-    smooth=False,
-    window_size=7,
-    poly_order=3,
-    variable="om",
-):
-
-    """function iterates through the dictionary created by load_cclv2 and locates peaks for each scan
-    args:   scan - a single scan,
-
-            int_threshold - fraction of threshold_intensity/max_intensity, must be positive num between 0 and 1
-                        i.e. will only detect peaks above 75% of max intensity
-
-            prominence - defines a drop of values that must be between two peaks, must be positive number
-                        i.e. if promimence is 20, it will detect two neigbouring peaks of 300 and 310 intesities,
-                        if none of the itermediate values are lower that 290
-
-            smooth  - if true, smooths data by savitzky golay filter, if false - no smoothing
-
-            window_size - window size for savgol filter, must be odd positive integer
-
-            poly_order =  order of the polynomial used in savgol filter, must be positive integer smaller than
-            window_size returns: dictionary with following structure:
-                                D{M34{  'num_of_peaks': 1,              #num of peaks
-                                        'peak_indexes': [20],           # index of peaks in omega array
-                                        'peak_heights': [90.],          # height of the peaks (if data vere smoothed
-                                                                        its the heigh of the peaks in smoothed data)
-    """
-    if not 0 <= int_threshold <= 1:
-        int_threshold = 0.8
-        print(
-            "Invalid value for int_threshold, select value between 0 and 1, new value set to:",
-            int_threshold,
-        )
-
-    if not isinstance(window_size, int) or (window_size % 2) == 0 or window_size <= 1:
-        window_size = 7
-        print(
-            "Invalid value for window_size, select positive odd integer, new value set to!:",
-            window_size,
-        )
-
-    if not isinstance(poly_order, int) or window_size < poly_order:
-        poly_order = 3
-        print(
-            "Invalid value for poly_order, select positive integer smaller than window_size, new value set to:",
-            poly_order,
-        )
-
-    if not isinstance(prominence, (int, float)) and prominence < 0:
-        prominence = 50
-        print("Invalid value for prominence, select positive number, new value set to:", prominence)
-
-    omega = scan[variable]
-    counts = np.array(scan["Counts"])
-    if smooth:
-        itp = interp1d(omega, counts, kind="linear")
-        absintensity = [abs(number) for number in counts]
-        lowest_intensity = min(absintensity)
-        counts[counts < 0] = lowest_intensity
-        smooth_peaks = savgol_filter(itp(omega), window_size, poly_order)
-
-    else:
-        smooth_peaks = counts
-
-    peaks, properties = sc.signal.find_peaks(
-        smooth_peaks, height=int_threshold * max(smooth_peaks), prominence=prominence
-    )
-    scan["num_of_peaks"] = len(peaks)
-    scan["peak_indexes"] = peaks
-    scan["peak_heights"] = properties["peak_heights"]
-    scan["smooth_peaks"] = smooth_peaks  # smoothed curve

pyzebra/ccl_io.py

@@ -1,5 +1,6 @@
 import os
 import re
+from ast import literal_eval
 from collections import defaultdict
 
 import numpy as np
@@ -55,41 +56,29 @@ META_VARS_FLOAT = (
     "s2hr",
     "s2hl",
 )
-META_UB_MATRIX = ("ub1j", "ub2j", "ub3j")
 
-CCL_FIRST_LINE = (
-    ("scan_number", int),
-    ("h_index", float),
-    ("k_index", float),
-    ("l_index", float),
-)
+META_UB_MATRIX = ("ub1j", "ub2j", "ub3j", "UB")
+
+CCL_FIRST_LINE = (("idx", int), ("h", float), ("k", float), ("l", float))
 
 CCL_ANGLES = {
-    "bi": (
-        ("twotheta_angle", float),
-        ("omega_angle", float),
-        ("chi_angle", float),
-        ("phi_angle", float),
-    ),
-    "nb": (
-        ("gamma_angle", float),
-        ("omega_angle", float),
-        ("nu_angle", float),
-        ("unkwn_angle", float),
-    ),
+    "bi": (("twotheta", float), ("omega", float), ("chi", float), ("phi", float)),
+    "nb": (("gamma", float), ("omega", float), ("nu", float), ("skip_angle", float)),
 }
 
 CCL_SECOND_LINE = (
     ("n_points", int),
     ("angle_step", float),
     ("monitor", float),
-    ("temperature", float),
-    ("mag_field", float),
+    ("temp", float),
+    ("mf", float),
     ("date", str),
     ("time", str),
-    ("scan_type", str),
+    ("scan_motor", str),
 )
 
+EXPORT_TARGETS = {"fullprof": (".comm", ".incomm"), "jana": (".col", ".incol")}
+
 
 def load_1D(filepath):
     """
@@ -105,67 +94,122 @@ def load_1D(filepath):
     """
     with open(filepath, "r") as infile:
         _, ext = os.path.splitext(filepath)
-        det_variables = parse_1D(infile, data_type=ext)
+        dataset = parse_1D(infile, data_type=ext)
 
-    return det_variables
+    return dataset
 
 
 def parse_1D(fileobj, data_type):
+    metadata = {"data_type": data_type}
+
     # read metadata
-    metadata = {}
     for line in fileobj:
         if "=" in line:
-            variable, value = line.split("=")
+            variable, value = line.split("=", 1)
             variable = variable.strip()
-            if variable in META_VARS_FLOAT:
-                metadata[variable] = float(value)
-            elif variable in META_VARS_STR:
-                metadata[variable] = str(value)[:-1].strip()
-            elif variable in META_UB_MATRIX:
-                metadata[variable] = re.findall(r"[-+]?\d*\.\d+|\d+", str(value))
+            value = value.strip()
+
+            try:
+                if variable in META_VARS_STR:
+                    metadata[variable] = value
+
+                elif variable in META_VARS_FLOAT:
+                    if variable == "2-theta":  # fix that angle name not to be an expression
+                        variable = "twotheta"
+                    if variable in ("a", "b", "c", "alpha", "beta", "gamma"):
+                        variable += "_cell"
+                    metadata[variable] = float(value)
+
+                elif variable in META_UB_MATRIX:
+                    if variable == "UB":
+                        metadata["ub"] = np.array(literal_eval(value)).reshape(3, 3)
+                    else:
+                        if "ub" not in metadata:
+                            metadata["ub"] = np.zeros((3, 3))
+                        row = int(variable[-2]) - 1
+                        metadata["ub"][row, :] = list(map(float, value.split()))
+
+            except Exception:
+                print(f"Error reading {variable} with value '{value}'")
+                metadata[variable] = 0
 
         if "#data" in line:
             # this is the end of metadata and the start of data section
             break
 
+    # handle older files that don't contain "zebra_mode" metadata
+    if "zebra_mode" not in metadata:
+        metadata["zebra_mode"] = "nb"
+
     # read data
-    scan = {}
+    dataset = []
     if data_type == ".ccl":
-        ccl_first_line = (*CCL_FIRST_LINE, *CCL_ANGLES[metadata["zebra_mode"]])
+        ccl_first_line = CCL_FIRST_LINE + CCL_ANGLES[metadata["zebra_mode"]]
         ccl_second_line = CCL_SECOND_LINE
 
         for line in fileobj:
-            d = {}
+            # skip empty/whitespace lines before start of any scan
+            if not line or line.isspace():
+                continue
+
+            scan = {}
+            scan["export"] = True
 
             # first line
             for param, (param_name, param_type) in zip(line.split(), ccl_first_line):
-                d[param_name] = param_type(param)
+                scan[param_name] = param_type(param)
 
             # second line
             next_line = next(fileobj)
             for param, (param_name, param_type) in zip(next_line.split(), ccl_second_line):
-                d[param_name] = param_type(param)
+                scan[param_name] = param_type(param)
 
-            d["om"] = np.linspace(
-                d["omega_angle"] - (d["n_points"] / 2) * d["angle_step"],
-                d["omega_angle"] + (d["n_points"] / 2) * d["angle_step"],
-                d["n_points"],
+            if "scan_motor" not in scan:
+                scan["scan_motor"] = "om"
+
+            if scan["scan_motor"] != "om":
+                raise Exception("Unsupported variable name in ccl file.")
+
+            # "om" -> "omega"
+            scan["scan_motor"] = "omega"
+            scan["scan_motors"] = ["omega", ]
+            # overwrite metadata, because it only refers to the scan center
+            half_dist = (scan["n_points"] - 1) / 2 * scan["angle_step"]
+            scan["omega"] = np.linspace(
+                scan["omega"] - half_dist, scan["omega"] + half_dist, scan["n_points"]
             )
 
             # subsequent lines with counts
             counts = []
-            while len(counts) < d["n_points"]:
-                counts.extend(map(int, next(fileobj).split()))
-            d["Counts"] = counts
+            while len(counts) < scan["n_points"]:
+                counts.extend(map(float, next(fileobj).split()))
+            scan["counts"] = np.array(counts)
+            scan["counts_err"] = np.sqrt(np.maximum(scan["counts"], 1))
 
-            scan[d["scan_number"]] = d
+            if scan["h"].is_integer() and scan["k"].is_integer() and scan["l"].is_integer():
+                scan["h"], scan["k"], scan["l"] = map(int, (scan["h"], scan["k"], scan["l"]))
+
+            dataset.append({**metadata, **scan})
 
     elif data_type == ".dat":
-        # skip the first 2 rows, the third row contans the column names
-        next(fileobj)
-        next(fileobj)
-        col_names = next(fileobj).split()
-        data_cols = defaultdict(list)
+        # TODO: this might need to be adapted in the future, when "gamma" will be added to dat files
+        if metadata["zebra_mode"] == "nb":
+            metadata["gamma"] = metadata["twotheta"]
+
+        scan = defaultdict(list)
+        scan["export"] = True
+
+        match = re.search("Scanning Variables: (.*), Steps: (.*)", next(fileobj))
+        motors = [motor.lower() for motor in match.group(1).split(", ")]
+        steps = [float(step) for step in match.group(2).split()]
+
+        match = re.search("(.*) Points, Mode: (.*), Preset (.*)", next(fileobj))
+        if match.group(2) != "Monitor":
+            raise Exception("Unknown mode in dat file.")
+        scan["n_points"] = int(match.group(1))
+        scan["monitor"] = float(match.group(3))
+
+        col_names = list(map(str.lower, next(fileobj).split()))
 
         for line in fileobj:
             if "END-OF-DATA" in line:
@@ -173,98 +217,193 @@ def parse_1D(fileobj, data_type):
                 break
 
             for name, val in zip(col_names, line.split()):
-                data_cols[name].append(float(val))
+                scan[name].append(float(val))
 
-        try:
-            data_cols["h_index"] = float(metadata["title"].split()[-3])
-            data_cols["k_index"] = float(metadata["title"].split()[-2])
-            data_cols["l_index"] = float(metadata["title"].split()[-1])
-        except (ValueError, IndexError):
-            print("seems hkl is not in title")
+        for name in col_names:
+            scan[name] = np.array(scan[name])
 
-        data_cols["temperature"] = metadata["temp"]
-        try:
-            data_cols["mag_field"] = metadata["mf"]
-        except KeyError:
-            print("Mag_field not present in dat file")
+        scan["counts_err"] = np.sqrt(np.maximum(scan["counts"], 1))
 
-        data_cols["omega_angle"] = metadata["omega"]
-        data_cols["n_points"] = len(data_cols["om"])
-        data_cols["monitor"] = data_cols["Monitor1"][0]
-        data_cols["twotheta_angle"] = metadata["2-theta"]
-        data_cols["chi_angle"] = metadata["chi"]
-        data_cols["phi_angle"] = metadata["phi"]
-        data_cols["nu_angle"] = metadata["nu"]
+        scan["scan_motors"] = []
+        for motor, step in zip(motors, steps):
+            if step == 0:
+                # it's not a scan motor, so keep only the median value
+                scan[motor] = np.median(scan[motor])
+            else:
+                scan["scan_motors"].append(motor)
 
-        scan[1] = dict(data_cols)
+        # "om" -> "omega"
+        if "om" in scan["scan_motors"]:
+            scan["scan_motors"][scan["scan_motors"].index("om")] = "omega"
+            scan["omega"] = scan["om"]
+            del scan["om"]
+
+        # "tt" -> "temp"
+        if "tt" in scan["scan_motors"]:
+            scan["scan_motors"][scan["scan_motors"].index("tt")] = "temp"
+            scan["temp"] = scan["tt"]
+            del scan["tt"]
+
+        # "mf" stays "mf"
+        # "phi" stays "phi"
+
+        scan["scan_motor"] = scan["scan_motors"][0]
+
+        if "h" not in scan:
+            scan["h"] = scan["k"] = scan["l"] = float("nan")
+
+        for param in ("mf", "temp"):
+            if param not in metadata:
+                scan[param] = 0
+
+        scan["idx"] = 1
+
+        dataset.append({**metadata, **scan})
 
     else:
         print("Unknown file extention")
 
-    # utility information
-    if all(
-        s["h_index"].is_integer() and s["k_index"].is_integer() and s["l_index"].is_integer()
-        for s in scan.values()
-    ):
-        metadata["indices"] = "hkl"
-    else:
-        metadata["indices"] = "real"
-
-    metadata["data_type"] = data_type
-    metadata["area_method"] = "fit"
-
-    return {"meta": metadata, "scan": scan}
+    return dataset
 
 
-def export_comm(data, path, lorentz=False):
-    """exports data in the *.comm format
-    :param lorentz: perform Lorentz correction
-    :param path: path to file + name
-    :arg data - data to export, is dict after peak fitting
+def export_1D(dataset, path, export_target, hkl_precision=2):
+    """Exports data in the .comm/.incomm format for fullprof or .col/.incol format for jana.
 
+    Scans with integer/real hkl values are saved in .comm/.incomm or .col/.incol files
+    correspondingly. If no scans are present for a particular output format, that file won't be
+    created.
     """
-    zebra_mode = data["meta"]["zebra_mode"]
-    if data["meta"]["indices"] == "hkl":
-        extension = ".comm"
-        padding = [6, 4]
-    elif data["meta"]["indices"] == "real":
-        extension = ".incomm"
-        padding = [4, 6]
+    if export_target not in EXPORT_TARGETS:
+        raise ValueError(f"Unknown export target: {export_target}.")
 
-    with open(str(path + extension), "w") as out_file:
-        for key, scan in data["scan"].items():
-            if "fit" not in scan:
-                print("Scan skipped - no fit value for:", key)
-                continue
+    zebra_mode = dataset[0]["zebra_mode"]
+    exts = EXPORT_TARGETS[export_target]
+    file_content = {ext: [] for ext in exts}
 
-            scan_str = f"{key:>{padding[0]}}"
-            h_str = f'{int(scan["h_index"]):{padding[1]}}'
-            k_str = f'{int(scan["k_index"]):{padding[1]}}'
-            l_str = f'{int(scan["l_index"]):{padding[1]}}'
+    for scan in dataset:
+        if "fit" not in scan:
+            continue
 
-            if data["meta"]["area_method"] == "fit":
-                area = scan["fit"]["fit_area"].n
-                sigma_str = f'{scan["fit"]["fit_area"].s:>10.2f}'
-            elif data["meta"]["area_method"] == "integ":
-                area = scan["fit"]["int_area"].n
-                sigma_str = f'{scan["fit"]["int_area"].s:>10.2f}'
+        idx_str = f"{scan['idx']:6}"
 
-            # apply lorentz correction to area
-            if lorentz:
-                if zebra_mode == "bi":
-                    twotheta_angle = np.deg2rad(scan["twotheta_angle"])
-                    corr_factor = np.sin(twotheta_angle)
-                elif zebra_mode == "nb":
-                    gamma_angle = np.deg2rad(scan["gamma_angle"])
-                    nu_angle = np.deg2rad(scan["nu_angle"])
-                    corr_factor = np.sin(gamma_angle) * np.cos(nu_angle)
+        h, k, l = scan["h"], scan["k"], scan["l"]
+        hkl_are_integers = isinstance(h, int)  # if True, other indices are of type 'int' too
+        if hkl_are_integers:
+            hkl_str = f"{h:4}{k:4}{l:4}"
+        else:
+            hkl_str = f"{h:8.{hkl_precision}f}{k:8.{hkl_precision}f}{l:8.{hkl_precision}f}"
 
-                area = np.abs(area * corr_factor)
+        area_n, area_s = scan["area"]
+        area_str = f"{area_n:10.2f}{area_s:10.2f}"
 
-            area_str = f"{area:>10.2f}"
+        ang_str = ""
+        for angle, _ in CCL_ANGLES[zebra_mode]:
+            if angle == scan["scan_motor"]:
+                angle_center = (np.min(scan[angle]) + np.max(scan[angle])) / 2
+            else:
+                angle_center = scan[angle]
 
-            ang_str = ""
-            for angle, _ in CCL_ANGLES[zebra_mode]:
-                ang_str = ang_str + f"{scan[angle]:8}"
+            if angle == "twotheta" and export_target == "jana":
+                angle_center /= 2
 
-            out_file.write(scan_str + h_str + k_str + l_str + area_str + sigma_str + ang_str + "\n")
+            ang_str = ang_str + f"{angle_center:8g}"
+
+        if export_target == "jana":
+            ang_str = ang_str + f"{scan['temp']:8}" + f"{scan['monitor']:8}"
+
+        ref = file_content[exts[0]] if hkl_are_integers else file_content[exts[1]]
+        ref.append(idx_str + hkl_str + area_str + ang_str + "\n")
+
+    for ext, content in file_content.items():
+        if content:
+            with open(path + ext, "w") as out_file:
+                out_file.writelines(content)
+
+
+def export_ccl_compare(dataset1, dataset2, path, export_target, hkl_precision=2):
+    """Exports compare data in the .comm/.incomm format for fullprof or .col/.incol format for jana.
+
+    Scans with integer/real hkl values are saved in .comm/.incomm or .col/.incol files
+    correspondingly. If no scans are present for a particular output format, that file won't be
+    created.
+    """
+    if export_target not in EXPORT_TARGETS:
+        raise ValueError(f"Unknown export target: {export_target}.")
+
+    zebra_mode = dataset1[0]["zebra_mode"]
+    exts = EXPORT_TARGETS[export_target]
+    file_content = {ext: [] for ext in exts}
+
+    for scan1, scan2 in zip(dataset1, dataset2):
+        if "fit" not in scan1:
+            continue
+
+        idx_str = f"{scan1['idx']:6}"
+
+        h, k, l = scan1["h"], scan1["k"], scan1["l"]
+        hkl_are_integers = isinstance(h, int)  # if True, other indices are of type 'int' too
+        if hkl_are_integers:
+            hkl_str = f"{h:4}{k:4}{l:4}"
+        else:
+            hkl_str = f"{h:8.{hkl_precision}f}{k:8.{hkl_precision}f}{l:8.{hkl_precision}f}"
+
+        area_n1, area_s1 = scan1["area"]
+        area_n2, area_s2 = scan2["area"]
+        area_n = area_n1 - area_n2
+        area_s = np.sqrt(area_s1 ** 2 + area_s2 ** 2)
+        area_str = f"{area_n:10.2f}{area_s:10.2f}"
+
+        ang_str = ""
+        for angle, _ in CCL_ANGLES[zebra_mode]:
+            if angle == scan1["scan_motor"]:
+                angle_center = (np.min(scan1[angle]) + np.max(scan1[angle])) / 2
+            else:
+                angle_center = scan1[angle]
+
+            if angle == "twotheta" and export_target == "jana":
+                angle_center /= 2
+
+            ang_str = ang_str + f"{angle_center:8g}"
+
+        if export_target == "jana":
+            ang_str = ang_str + f"{scan1['temp']:8}" + f"{scan1['monitor']:8}"
+
+        ref = file_content[exts[0]] if hkl_are_integers else file_content[exts[1]]
+        ref.append(idx_str + hkl_str + area_str + ang_str + "\n")
+
+    for ext, content in file_content.items():
+        if content:
+            with open(path + ext, "w") as out_file:
+                out_file.writelines(content)
+
+
+def export_param_study(dataset, param_data, path):
+    file_content = []
+    for scan, param in zip(dataset, param_data):
+        if "fit" not in scan:
+            continue
+
+        if not file_content:
+            title_str = f"{'param':12}"
+            for fit_param_name in scan["fit"].params:
+                title_str = title_str + f"{fit_param_name:20}" + f"{'std_' + fit_param_name:20}"
+            title_str = title_str + "file"
+            file_content.append(title_str + "\n")
+
+        param_str = f"{param:<12.2f}"
+
+        fit_str = ""
+        for fit_param in scan["fit"].params.values():
+            fit_param_val = fit_param.value
+            fit_param_std = fit_param.stderr
+            if fit_param_std is None:
+                fit_param_std = 0
+            fit_str = fit_str + f"{fit_param_val:<20.2f}" + f"{fit_param_std:<20.2f}"
+
+        _, fname_str = os.path.split(scan["original_filename"])
+
+        file_content.append(param_str + fit_str + fname_str + "\n")
+
+    if file_content:
+        with open(path, "w") as out_file:
+            out_file.writelines(file_content)

pyzebra/ccl_process.py

@@ -0,0 +1,338 @@
+import os
+
+import numpy as np
+from lmfit.models import GaussianModel, LinearModel, PseudoVoigtModel, VoigtModel
+from scipy.integrate import simpson, trapezoid
+
+from .ccl_io import CCL_ANGLES
+
+PARAM_PRECISIONS = {
+    "twotheta": 0.1,
+    "chi": 0.1,
+    "nu": 0.1,
+    "phi": 0.05,
+    "omega": 0.05,
+    "gamma": 0.05,
+    "temp": 1,
+    "mf": 0.001,
+    "ub": 0.01,
+}
+
+MAX_RANGE_GAP = {
+    "omega": 0.5,
+}
+
+MOTOR_POS_PRECISION = 0.01
+
+AREA_METHODS = ("fit_area", "int_area")
+
+
+def normalize_dataset(dataset, monitor=100_000):
+    for scan in dataset:
+        monitor_ratio = monitor / scan["monitor"]
+        scan["counts"] *= monitor_ratio
+        scan["counts_err"] *= monitor_ratio
+        scan["monitor"] = monitor
+
+
+def merge_duplicates(dataset):
+    merged = np.zeros(len(dataset), dtype=np.bool)
+    for ind_into, scan_into in enumerate(dataset):
+        for ind_from, scan_from in enumerate(dataset[ind_into + 1 :], start=ind_into + 1):
+            if _parameters_match(scan_into, scan_from) and not merged[ind_from]:
+                merge_scans(scan_into, scan_from)
+                merged[ind_from] = True
+
+
+def _parameters_match(scan1, scan2):
+    zebra_mode = scan1["zebra_mode"]
+    if zebra_mode != scan2["zebra_mode"]:
+        return False
+
+    for param in ("ub", *(vars[0] for vars in CCL_ANGLES[zebra_mode])):
+        if param.startswith("skip"):
+            # ignore skip parameters, like the last angle in 'nb' zebra mode
+            continue
+
+        if param == scan1["scan_motor"] == scan2["scan_motor"]:
+            # check if ranges of variable parameter overlap
+            r1_start, r1_end = scan1[param][0], scan1[param][-1]
+            r2_start, r2_end = scan2[param][0], scan2[param][-1]
+            # support reversed ranges
+            if r1_start > r1_end:
+                r1_start, r1_end = r1_end, r1_start
+            if r2_start > r2_end:
+                r2_start, r2_end = r2_end, r2_start
+            # maximum gap between ranges of the scanning parameter (default 0)
+            max_range_gap = MAX_RANGE_GAP.get(param, 0)
+            if max(r1_start - r2_end, r2_start - r1_end) > max_range_gap:
+                return False
+
+        elif (
+            np.max(np.abs(np.median(scan1[param]) - np.median(scan2[param])))
+            > PARAM_PRECISIONS[param]
+        ):
+            return False
+
+    return True
+
+
+def merge_datasets(dataset_into, dataset_from):
+    scan_motors_into = dataset_into[0]["scan_motors"]
+    scan_motors_from = dataset_from[0]["scan_motors"]
+    if scan_motors_into != scan_motors_from:
+        print(f"Scan motors mismatch between datasets: {scan_motors_into} vs {scan_motors_from}")
+        return
+
+    merged = np.zeros(len(dataset_from), dtype=np.bool)
+    for scan_into in dataset_into:
+        for ind, scan_from in enumerate(dataset_from):
+            if _parameters_match(scan_into, scan_from) and not merged[ind]:
+                if scan_into["counts"].ndim == 3:
+                    merge_h5_scans(scan_into, scan_from)
+                else:  # scan_into["counts"].ndim == 1
+                    merge_scans(scan_into, scan_from)
+                merged[ind] = True
+
+    for scan_from in dataset_from:
+        dataset_into.append(scan_from)
+
+
+def merge_scans(scan_into, scan_from):
+    if "init_scan" not in scan_into:
+        scan_into["init_scan"] = scan_into.copy()
+
+    if "merged_scans" not in scan_into:
+        scan_into["merged_scans"] = []
+
+    if scan_from in scan_into["merged_scans"]:
+        return
+
+    scan_into["merged_scans"].append(scan_from)
+
+    scan_motor = scan_into["scan_motor"]  # the same as scan_from["scan_motor"]
+
+    pos_all = np.array([])
+    val_all = np.array([])
+    err_all = np.array([])
+    for scan in [scan_into["init_scan"], *scan_into["merged_scans"]]:
+        pos_all = np.append(pos_all, scan[scan_motor])
+        val_all = np.append(val_all, scan["counts"])
+        err_all = np.append(err_all, scan["counts_err"] ** 2)
+
+    sort_index = np.argsort(pos_all)
+    pos_all = pos_all[sort_index]
+    val_all = val_all[sort_index]
+    err_all = err_all[sort_index]
+
+    pos_tmp = pos_all[:1]
+    val_tmp = val_all[:1]
+    err_tmp = err_all[:1]
+    num_tmp = np.array([1])
+    for pos, val, err in zip(pos_all[1:], val_all[1:], err_all[1:]):
+        if pos - pos_tmp[-1] < MOTOR_POS_PRECISION:
+            # the repeated motor position
+            val_tmp[-1] += val
+            err_tmp[-1] += err
+            num_tmp[-1] += 1
+        else:
+            # a new motor position
+            pos_tmp = np.append(pos_tmp, pos)
+            val_tmp = np.append(val_tmp, val)
+            err_tmp = np.append(err_tmp, err)
+            num_tmp = np.append(num_tmp, 1)
+
+    scan_into[scan_motor] = pos_tmp
+    scan_into["counts"] = val_tmp / num_tmp
+    scan_into["counts_err"] = np.sqrt(err_tmp) / num_tmp
+
+    scan_from["export"] = False
+
+    fname1 = os.path.basename(scan_into["original_filename"])
+    fname2 = os.path.basename(scan_from["original_filename"])
+    print(f'Merging scans: {scan_into["idx"]} ({fname1}) <-- {scan_from["idx"]} ({fname2})')
+
+
+def merge_h5_scans(scan_into, scan_from):
+    if "init_scan" not in scan_into:
+        scan_into["init_scan"] = scan_into.copy()
+
+    if "merged_scans" not in scan_into:
+        scan_into["merged_scans"] = []
+
+    for scan in scan_into["merged_scans"]:
+        if scan_from is scan:
+            print("Already merged scan")
+            return
+
+    scan_into["merged_scans"].append(scan_from)
+
+    scan_motor = scan_into["scan_motor"]  # the same as scan_from["scan_motor"]
+
+    pos_all = [scan_into["init_scan"][scan_motor]]
+    val_all = [scan_into["init_scan"]["counts"]]
+    err_all = [scan_into["init_scan"]["counts_err"] ** 2]
+    for scan in scan_into["merged_scans"]:
+        pos_all.append(scan[scan_motor])
+        val_all.append(scan["counts"])
+        err_all.append(scan["counts_err"] ** 2)
+    pos_all = np.concatenate(pos_all)
+    val_all = np.concatenate(val_all)
+    err_all = np.concatenate(err_all)
+
+    sort_index = np.argsort(pos_all)
+    pos_all = pos_all[sort_index]
+    val_all = val_all[sort_index]
+    err_all = err_all[sort_index]
+
+    pos_tmp = [pos_all[0]]
+    val_tmp = [val_all[:1]]
+    err_tmp = [err_all[:1]]
+    num_tmp = [1]
+    for pos, val, err in zip(pos_all[1:], val_all[1:], err_all[1:]):
+        if pos - pos_tmp[-1] < MOTOR_POS_PRECISION:
+            # the repeated motor position
+            val_tmp[-1] += val
+            err_tmp[-1] += err
+            num_tmp[-1] += 1
+        else:
+            # a new motor position
+            pos_tmp.append(pos)
+            val_tmp.append(val[None, :])
+            err_tmp.append(err[None, :])
+            num_tmp.append(1)
+    pos_tmp = np.array(pos_tmp)
+    val_tmp = np.concatenate(val_tmp)
+    err_tmp = np.concatenate(err_tmp)
+    num_tmp = np.array(num_tmp)
+
+    scan_into[scan_motor] = pos_tmp
+    scan_into["counts"] = val_tmp / num_tmp[:, None, None]
+    scan_into["counts_err"] = np.sqrt(err_tmp) / num_tmp[:, None, None]
+
+    scan_from["export"] = False
+
+    fname1 = os.path.basename(scan_into["original_filename"])
+    fname2 = os.path.basename(scan_from["original_filename"])
+    print(f'Merging scans: {scan_into["idx"]} ({fname1}) <-- {scan_from["idx"]} ({fname2})')
+
+
+def restore_scan(scan):
+    if "merged_scans" in scan:
+        for merged_scan in scan["merged_scans"]:
+            merged_scan["export"] = True
+
+    if "init_scan" in scan:
+        tmp = scan["init_scan"]
+        scan.clear()
+        scan.update(tmp)
+        # force scan export to True, otherwise in the sequence of incorrectly merged scans
+        # a <- b <- c the scan b will be restored with scan["export"] = False if restoring executed
+        # in the same order, i.e. restore a -> restore b
+        scan["export"] = True
+
+
+def fit_scan(scan, model_dict, fit_from=None, fit_to=None):
+    if fit_from is None:
+        fit_from = -np.inf
+    if fit_to is None:
+        fit_to = np.inf
+
+    y_fit = scan["counts"]
+    y_err = scan["counts_err"]
+    x_fit = scan[scan["scan_motor"]]
+
+    # apply fitting range
+    fit_ind = (fit_from <= x_fit) & (x_fit <= fit_to)
+    if not np.any(fit_ind):
+        print(f"No data in fit range for scan {scan['idx']}")
+        return
+
+    y_fit = y_fit[fit_ind]
+    y_err = y_err[fit_ind]
+    x_fit = x_fit[fit_ind]
+
+    model = None
+    for model_index, (model_name, model_param) in enumerate(model_dict.items()):
+        model_name, _ = model_name.split("-")
+        prefix = f"f{model_index}_"
+
+        if model_name == "linear":
+            _model = LinearModel(prefix=prefix)
+        elif model_name == "gaussian":
+            _model = GaussianModel(prefix=prefix)
+        elif model_name == "voigt":
+            _model = VoigtModel(prefix=prefix)
+        elif model_name == "pvoigt":
+            _model = PseudoVoigtModel(prefix=prefix)
+        else:
+            raise ValueError(f"Unknown model name: '{model_name}'")
+
+        _init_guess = _model.guess(y_fit, x=x_fit)
+
+        for param_index, param_name in enumerate(model_param["param"]):
+            param_hints = {}
+            for hint_name in ("value", "vary", "min", "max"):
+                tmp = model_param[hint_name][param_index]
+                if tmp is None:
+                    param_hints[hint_name] = getattr(_init_guess[prefix + param_name], hint_name)
+                else:
+                    param_hints[hint_name] = tmp
+
+            if "center" in param_name:
+                if np.isneginf(param_hints["min"]):
+                    param_hints["min"] = np.min(x_fit)
+
+                if np.isposinf(param_hints["max"]):
+                    param_hints["max"] = np.max(x_fit)
+
+            if "sigma" in param_name:
+                if np.isposinf(param_hints["max"]):
+                    param_hints["max"] = np.max(x_fit) - np.min(x_fit)
+
+            _model.set_param_hint(param_name, **param_hints)
+
+        if model is None:
+            model = _model
+        else:
+            model += _model
+
+    scan["fit"] = model.fit(y_fit, x=x_fit, weights=1 / y_err)
+
+
+def get_area(scan, area_method, lorentz):
+    if "fit" not in scan:
+        return
+
+    if area_method not in AREA_METHODS:
+        raise ValueError(f"Unknown area method: {area_method}.")
+
+    if area_method == "fit_area":
+        area_v = 0
+        area_s = 0
+        for name, param in scan["fit"].params.items():
+            if "amplitude" in name:
+                area_v += np.nan if param.value is None else param.value
+                area_s += np.nan if param.stderr is None else param.stderr
+
+    else:  # area_method == "int_area"
+        y_val = scan["counts"]
+        x_val = scan[scan["scan_motor"]]
+        y_bkg = scan["fit"].eval_components(x=x_val)["f0_"]
+        area_v = simpson(y_val, x=x_val) - trapezoid(y_bkg, x=x_val)
+        area_s = np.sqrt(area_v)
+
+    if lorentz:
+        # lorentz correction to area
+        if scan["zebra_mode"] == "bi":
+            twotheta = np.deg2rad(scan["twotheta"])
+            corr_factor = np.sin(twotheta)
+        else:  # zebra_mode == "nb":
+            gamma = np.deg2rad(scan["gamma"])
+            nu = np.deg2rad(scan["nu"])
+            corr_factor = np.sin(gamma) * np.cos(nu)
+
+        area_v = np.abs(area_v * corr_factor)
+        area_s = np.abs(area_s * corr_factor)
+
+    scan["area"] = (area_v, area_s)

pyzebra/fit2.py

@@ -1,228 +0,0 @@
-import numpy as np
-import uncertainties as u
-from lmfit import Model, Parameters
-from scipy.integrate import simps
-
-
-def bin_data(array, binsize):
-    if isinstance(binsize, int) and 0 < binsize < len(array):
-        return [
-            np.mean(array[binsize * i : binsize * i + binsize])
-            for i in range(int(np.ceil(len(array) / binsize)))
-        ]
-    else:
-        print("Binsize need to be positive integer smaller than lenght of array")
-        return array
-
-
-def find_nearest(array, value):
-    # find nearest value and return index
-    array = np.asarray(array)
-    idx = (np.abs(array - value)).argmin()
-    return idx
-
-
-def create_uncertanities(y, y_err):
-    # create array with uncertanities for error propagation
-    combined = np.array([])
-    for i in range(len(y)):
-        part = u.ufloat(y[i], y_err[i])
-        combined = np.append(combined, part)
-    return combined
-
-
-def fitccl(
-    scan,
-    guess,
-    vary,
-    constraints_min,
-    constraints_max,
-    numfit_min=None,
-    numfit_max=None,
-    binning=None,
-):
-    """Made for fitting of ccl date where 1 peak is expected. Allows for combination of gaussian and linear model combination
-    :param scan: scan in the data dict (i.e. M123)
-    :param guess: initial guess for the fitting, if none, some values are added automatically in order (see below)
-    :param vary: True if parameter can vary during fitting, False if it to be fixed
-    :param numfit_min: minimal value on x axis for numerical integration - if none is centre of gaussian minus 3 sigma
-    :param numfit_max: maximal value on x axis for numerical integration - if none is centre of gaussian plus 3 sigma
-    :param constraints_min: min constranits value for fit
-    :param constraints_max: max constranits value for fit
-    :param binning : binning of the data
-    :return data dict with additional values
-    order for guess, vary, constraints_min, constraints_max:
-    [Gaussian centre, Gaussian sigma, Gaussian amplitude, background slope, background intercept]
-    examples:
-    guess = [None, None, 100,  0, None]
-    vary = [True, True, True, True,  True]
-    constraints_min = [23, None, 50, 0, 0]
-    constraints_min = [80, None, 1000, 0, 100]
-    """
-    if "peak_indexes" not in scan:
-        scan["peak_indexes"] = []
-
-    if len(scan["peak_indexes"]) > 1:
-        # return in case of more than 1 peaks
-        return
-
-    if binning is None or binning == 0 or binning == 1:
-        x = list(scan["om"])
-        y = list(scan["Counts"])
-        y_err = list(np.sqrt(y)) if scan.get("sigma", None) is None else list(scan["sigma"])
-        if not scan["peak_indexes"]:
-            centre = np.mean(x)
-        else:
-            centre = x[int(scan["peak_indexes"])]
-    else:
-        x = list(scan["om"])
-        if not scan["peak_indexes"]:
-            centre = np.mean(x)
-        else:
-            centre = x[int(scan["peak_indexes"])]
-        x = bin_data(x, binning)
-        y = list(scan["Counts"])
-        y_err = list(np.sqrt(y)) if scan.get("sigma", None) is None else list(scan["sigma"])
-        combined = bin_data(create_uncertanities(y, y_err), binning)
-        y = [combined[i].n for i in range(len(combined))]
-        y_err = [combined[i].s for i in range(len(combined))]
-
-    if len(scan["peak_indexes"]) == 0:
-        # Case for no peak, gaussian in centre, sigma as 20% of range
-        peak_index = find_nearest(x, np.mean(x))
-        guess[0] = centre if guess[0] is None else guess[0]
-        guess[1] = (x[-1] - x[0]) / 5 if guess[1] is None else guess[1]
-        guess[2] = 50 if guess[2] is None else guess[2]
-        guess[3] = 0 if guess[3] is None else guess[3]
-        guess[4] = np.mean(y) if guess[4] is None else guess[4]
-        constraints_min[2] = 0
-
-    elif len(scan["peak_indexes"]) == 1:
-        # case for one peak, takse into account users guesses
-        peak_height = scan["peak_heights"]
-        guess[0] = centre if guess[0] is None else guess[0]
-        guess[1] = 0.1 if guess[1] is None else guess[1]
-        guess[2] = float(peak_height / 10) if guess[2] is None else float(guess[2])
-        guess[3] = 0 if guess[3] is None else guess[3]
-        guess[4] = np.median(x) if guess[4] is None else guess[4]
-        constraints_min[0] = np.min(x) if constraints_min[0] is None else constraints_min[0]
-        constraints_max[0] = np.max(x) if constraints_max[0] is None else constraints_max[0]
-
-    def gaussian(x, g_cen, g_width, g_amp):
-        """1-d gaussian: gaussian(x, amp, cen, wid)"""
-        return (g_amp / (np.sqrt(2 * np.pi) * g_width)) * np.exp(
-            -((x - g_cen) ** 2) / (2 * g_width ** 2)
-        )
-
-    def background(x, slope, intercept):
-        """background"""
-        return slope * (x - centre) + intercept
-
-    mod = Model(gaussian) + Model(background)
-    params = Parameters()
-    params.add_many(
-        ("g_cen", guess[0], bool(vary[0]), np.min(x), np.max(x), None, None),
-        ("g_width", guess[1], bool(vary[1]), constraints_min[1], constraints_max[1], None, None),
-        ("g_amp", guess[2], bool(vary[2]), constraints_min[2], constraints_max[2], None, None),
-        ("slope", guess[3], bool(vary[3]), constraints_min[3], constraints_max[3], None, None),
-        ("intercept", guess[4], bool(vary[4]), constraints_min[4], constraints_max[4], None, None),
-    )
-    # the weighted fit
-    weights = [np.abs(1 / val) if val != 0 else 1 for val in y_err]
-    try:
-        result = mod.fit(y, params, weights=weights, x=x, calc_covar=True)
-    except ValueError:
-        print(f"Couldn't fit scan {scan['scan_number']}")
-        return
-
-    if result.params["g_amp"].stderr is None:
-        result.params["g_amp"].stderr = result.params["g_amp"].value
-    elif result.params["g_amp"].stderr > result.params["g_amp"].value:
-        result.params["g_amp"].stderr = result.params["g_amp"].value
-
-    # u.ufloat to work with uncertanities
-    fit_area = u.ufloat(result.params["g_amp"].value, result.params["g_amp"].stderr)
-    comps = result.eval_components()
-
-    if len(scan["peak_indexes"]) == 0:
-        # for case of no peak, there is no reason to integrate, therefore fit and int are equal
-        int_area = fit_area
-
-    elif len(scan["peak_indexes"]) == 1:
-        gauss_3sigmamin = find_nearest(
-            x, result.params["g_cen"].value - 3 * result.params["g_width"].value
-        )
-        gauss_3sigmamax = find_nearest(
-            x, result.params["g_cen"].value + 3 * result.params["g_width"].value
-        )
-        numfit_min = gauss_3sigmamin if numfit_min is None else find_nearest(x, numfit_min)
-        numfit_max = gauss_3sigmamax if numfit_max is None else find_nearest(x, numfit_max)
-
-        it = -1
-        while abs(numfit_max - numfit_min) < 3:
-            # in the case the peak is very thin and numerical integration would be on zero omega
-            # difference, finds closes values
-            it = it + 1
-            numfit_min = find_nearest(
-                x,
-                result.params["g_cen"].value - 3 * (1 + it / 10) * result.params["g_width"].value,
-            )
-            numfit_max = find_nearest(
-                x,
-                result.params["g_cen"].value + 3 * (1 + it / 10) * result.params["g_width"].value,
-            )
-
-        if x[numfit_min] < np.min(x):
-            # makes sure that the values supplied by user lay in the omega range
-            # can be ommited for users who know what they're doing
-            numfit_min = gauss_3sigmamin
-            print("Minimal integration value outside of x range")
-        elif x[numfit_min] >= x[numfit_max]:
-            numfit_min = gauss_3sigmamin
-            print("Minimal integration value higher than maximal")
-        else:
-            pass
-        if x[numfit_max] > np.max(x):
-            numfit_max = gauss_3sigmamax
-            print("Maximal integration value outside of x range")
-        elif x[numfit_max] <= x[numfit_min]:
-            numfit_max = gauss_3sigmamax
-            print("Maximal integration value lower than minimal")
-        else:
-            pass
-
-        count_errors = create_uncertanities(y, y_err)
-        # create error vector for numerical integration propagation
-        num_int_area = simps(count_errors[numfit_min:numfit_max], x[numfit_min:numfit_max])
-        slope_err = u.ufloat(result.params["slope"].value, result.params["slope"].stderr)
-        # pulls the nominal and error values from fit (slope)
-        intercept_err = u.ufloat(
-            result.params["intercept"].value, result.params["intercept"].stderr
-        )
-        # pulls the nominal and error values from fit (intercept)
-
-        background_errors = np.array([])
-        for j in range(len(x[numfit_min:numfit_max])):
-            # creates nominal and error vector for numerical integration of background
-            bg = slope_err * (x[j] - centre) + intercept_err
-            background_errors = np.append(background_errors, bg)
-
-        num_int_background = simps(background_errors, x[numfit_min:numfit_max])
-        int_area = num_int_area - num_int_background
-
-    d = {}
-    for pars in result.params:
-        d[str(pars)] = (result.params[str(pars)].value, result.params[str(pars)].vary)
-
-    print("Scan", scan["scan_number"])
-    print(result.fit_report())
-
-    d["ratio"] = (result.params["g_amp"].value - int_area.n) / result.params["g_amp"].value
-    d["int_area"] = int_area
-    d["fit_area"] = u.ufloat(result.params["g_amp"].value, result.params["g_amp"].stderr)
-    d["full_report"] = result.fit_report()
-    d["result"] = result
-    d["comps"] = comps
-    d["numfit"] = [numfit_min, numfit_max]
-    d["x_fit"] = x
-    scan["fit"] = d

pyzebra/fitvol3.py

@@ -1,167 +0,0 @@
-import numpy as np
-from lmfit import Model, Parameters
-from scipy.integrate import simps
-import matplotlib.pyplot as plt
-import uncertainties as u
-from lmfit.models import GaussianModel
-from lmfit.models import VoigtModel
-from lmfit.models import PseudoVoigtModel
-
-
-def bin_data(array, binsize):
-    if isinstance(binsize, int) and 0 < binsize < len(array):
-        return [
-            np.mean(array[binsize * i : binsize * i + binsize])
-            for i in range(int(np.ceil(len(array) / binsize)))
-        ]
-    else:
-        print("Binsize need to be positive integer smaller than lenght of array")
-        return array
-
-
-def create_uncertanities(y, y_err):
-    # create array with uncertanities for error propagation
-    combined = np.array([])
-    for i in range(len(y)):
-        part = u.ufloat(y[i], y_err[i])
-        combined = np.append(combined, part)
-    return combined
-
-
-def find_nearest(array, value):
-    # find nearest value and return index
-    array = np.asarray(array)
-    idx = (np.abs(array - value)).argmin()
-    return idx
-
-
-# predefined peak positions
-# peaks = [6.2,  8.1, 9.9, 11.5]
-peaks = [23.5, 24.5]
-# peaks = [24]
-def fitccl(scan, variable="om", peak_type="gauss", binning=None):
-
-    x = list(scan[variable])
-    y = list(scan["Counts"])
-    peak_centre = np.mean(x)
-    if binning is None or binning == 0 or binning == 1:
-        x = list(scan["om"])
-        y = list(scan["Counts"])
-        y_err = list(np.sqrt(y)) if scan.get("sigma", None) is None else list(scan["sigma"])
-        print(scan["peak_indexes"])
-        if not scan["peak_indexes"]:
-            peak_centre = np.mean(x)
-        else:
-            centre = x[int(scan["peak_indexes"])]
-    else:
-        x = list(scan["om"])
-        if not scan["peak_indexes"]:
-            peak_centre = np.mean(x)
-        else:
-            peak_centre = x[int(scan["peak_indexes"])]
-        x = bin_data(x, binning)
-        y = list(scan["Counts"])
-        y_err = list(np.sqrt(y)) if scan.get("sigma", None) is None else list(scan["sigma"])
-        combined = bin_data(create_uncertanities(y, y_err), binning)
-        y = [combined[i].n for i in range(len(combined))]
-        y_err = [combined[i].s for i in range(len(combined))]
-
-    def background(x, slope, intercept):
-        """background"""
-        return slope * (x - peak_centre) + intercept
-
-    def gaussian(x, center, g_sigma, amplitude):
-        """1-d gaussian: gaussian(x, amp, cen, wid)"""
-        return (amplitude / (np.sqrt(2.0 * np.pi) * g_sigma)) * np.exp(
-            -((x - center) ** 2) / (2 * g_sigma ** 2)
-        )
-
-    def lorentzian(x, center, l_sigma, amplitude):
-        """1d lorentzian"""
-        return (amplitude / (1 + ((1 * x - center) / l_sigma) ** 2)) / (np.pi * l_sigma)
-
-    def pseudoVoigt1(x, center, g_sigma, amplitude, l_sigma, fraction):
-        """PseudoVoight peak with different widths of lorenzian and gaussian"""
-        return (1 - fraction) * gaussian(x, center, g_sigma, amplitude) + fraction * (
-            lorentzian(x, center, l_sigma, amplitude)
-        )
-
-    mod = Model(background)
-    params = Parameters()
-    params.add_many(
-        ("slope", 0, True, None, None, None, None), ("intercept", 0, False, None, None, None, None)
-    )
-    for i in range(len(peaks)):
-        if peak_type == "gauss":
-            mod = mod + GaussianModel(prefix="p%d_" % (i + 1))
-            params.add(str("p%d_" % (i + 1) + "amplitude"), 20, True, 0, None, None)
-            params.add(str("p%d_" % (i + 1) + "center"), peaks[i], True, None, None, None)
-            params.add(str("p%d_" % (i + 1) + "sigma"), 0.2, True, 0, 5, None)
-        elif peak_type == "voigt":
-            mod = mod + VoigtModel(prefix="p%d_" % (i + 1))
-            params.add(str("p%d_" % (i + 1) + "amplitude"), 20, True, 0, None, None)
-            params.add(str("p%d_" % (i + 1) + "center"), peaks[i], True, None, None, None)
-            params.add(str("p%d_" % (i + 1) + "sigma"), 0.2, True, 0, 3, None)
-            params.add(str("p%d_" % (i + 1) + "gamma"), 0.2, True, 0, 5, None)
-        elif peak_type == "pseudovoigt":
-            mod = mod + PseudoVoigtModel(prefix="p%d_" % (i + 1))
-            params.add(str("p%d_" % (i + 1) + "amplitude"), 20, True, 0, None, None)
-            params.add(str("p%d_" % (i + 1) + "center"), peaks[i], True, None, None, None)
-            params.add(str("p%d_" % (i + 1) + "sigma"), 0.2, True, 0, 5, None)
-            params.add(str("p%d_" % (i + 1) + "fraction"), 0.5, True, -5, 5, None)
-        elif peak_type == "pseudovoigt1":
-            mod = mod + Model(pseudoVoigt1, prefix="p%d_" % (i + 1))
-            params.add(str("p%d_" % (i + 1) + "amplitude"), 20, True, 0, None, None)
-            params.add(str("p%d_" % (i + 1) + "center"), peaks[i], True, None, None, None)
-            params.add(str("p%d_" % (i + 1) + "g_sigma"), 0.2, True, 0, 5, None)
-            params.add(str("p%d_" % (i + 1) + "l_sigma"), 0.2, True, 0, 5, None)
-            params.add(str("p%d_" % (i + 1) + "fraction"), 0.5, True, 0, 1, None)
-    # add parameters
-
-    result = mod.fit(
-        y, params, weights=[np.abs(1 / y_err[i]) for i in range(len(y_err))], x=x, calc_covar=True
-    )
-
-    comps = result.eval_components()
-
-    reportstring = list()
-    for keys in result.params:
-        if result.params[keys].value is not None:
-            str2 = np.around(result.params[keys].value, 3)
-        else:
-            str2 = 0
-        if result.params[keys].stderr is not None:
-            str3 = np.around(result.params[keys].stderr, 3)
-        else:
-            str3 = 0
-        reportstring.append("%s = %2.3f +/- %2.3f" % (keys, str2, str3))
-
-    reportstring = "\n".join(reportstring)
-
-    plt.figure(figsize=(20, 10))
-    plt.plot(x, result.best_fit, "k-", label="Best fit")
-
-    plt.plot(x, y, "b-", label="Original data")
-    plt.plot(x, comps["background"], "g--", label="Line component")
-    for i in range(len(peaks)):
-        plt.plot(
-            x,
-            comps[str("p%d_" % (i + 1))],
-            "r--",
-        )
-        plt.fill_between(x, comps[str("p%d_" % (i + 1))], alpha=0.4, label=str("p%d_" % (i + 1)))
-    plt.legend()
-    plt.text(
-        np.min(x),
-        np.max(y),
-        reportstring,
-        fontsize=9,
-        verticalalignment="top",
-    )
-    plt.title(str(peak_type))
-
-    plt.xlabel("Omega [deg]")
-    plt.ylabel("Counts [a.u.]")
-    plt.show()
-
-    print(result.fit_report())

pyzebra/h5.py

@@ -1,5 +1,10 @@
 import h5py
+import numpy as np
+from lmfit.models import Gaussian2dModel, GaussianModel
 
+META_MATRIX = ("UB", )
+META_CELL = ("cell", )
+META_STR = ("name", )
 
 def read_h5meta(filepath):
     """Open and parse content of a h5meta file.
@@ -23,49 +28,167 @@ def parse_h5meta(file):
         line = line.strip()
         if line.startswith("#begin "):
             section = line[len("#begin ") :]
-            content[section] = []
+            if section in ("detector parameters", "crystal"):
+                content[section] = {}
+            else:
+                content[section] = []
 
         elif line.startswith("#end"):
             section = None
 
         elif section:
-            content[section].append(line)
+            if section in ("detector parameters", "crystal"):
+                if "=" in line:
+                    variable, value = line.split("=", 1)
+                    variable = variable.strip()
+                    value = value.strip()
+
+                    if variable in META_STR:
+                        pass
+                    elif variable in META_CELL:
+                        value = np.array(value.split(",")[:6], dtype=np.float)
+                    elif variable in META_MATRIX:
+                        value = np.array(value.split(",")[:9], dtype=np.float).reshape(3, 3)
+                    else:  # default is a single float number
+                        value = float(value)
+                    content[section][variable] = value
+            else:
+                content[section].append(line)
 
     return content
 
 
-def read_detector_data(filepath):
+def read_detector_data(filepath, cami_meta=None):
     """Read detector data and angles from an h5 file.
 
     Args:
         filepath (str): File path of an h5 file.
 
     Returns:
-        ndarray: A 3D array of data, rot_angle, pol_angle, tilt_angle.
+        ndarray: A 3D array of data, omega, gamma, nu.
     """
     with h5py.File(filepath, "r") as h5f:
-        data = h5f["/entry1/area_detector2/data"][:]
+        counts = h5f["/entry1/area_detector2/data"][:].astype(np.float64)
 
-        # reshape data to a correct shape (2006 issue)
-        n, cols, rows = data.shape
-        data = data.reshape(n, rows, cols)
+        n, cols, rows = counts.shape
+        if "/entry1/experiment_identifier" in h5f:  # old format
+            # reshape images (counts) to a correct shape (2006 issue)
+            counts = counts.reshape(n, rows, cols)
+        else:
+            counts = counts.swapaxes(1, 2)
 
-        det_data = {"data": data}
+        scan = {"counts": counts, "counts_err": np.sqrt(np.maximum(counts, 1))}
+        scan["original_filename"] = filepath
+        scan["export"] = True
 
-        det_data["rot_angle"] = h5f["/entry1/area_detector2/rotation_angle"][:]  # om, sometimes ph
-        det_data["pol_angle"] = h5f["/entry1/ZEBRA/area_detector2/polar_angle"][:]  # gammad
-        det_data["tlt_angle"] = h5f["/entry1/ZEBRA/area_detector2/tilt_angle"][:]  # nud
-        det_data["ddist"] = h5f["/entry1/ZEBRA/area_detector2/distance"][:]
-        det_data["wave"] = h5f["/entry1/ZEBRA/monochromator/wavelength"][:]
-        det_data["chi_angle"] = h5f["/entry1/sample/chi"][:]  # ch
-        det_data["phi_angle"] = h5f["/entry1/sample/phi"][:]  # ph
-        det_data["UB"] = h5f["/entry1/sample/UB"][:].reshape(3, 3)
+        if "/entry1/zebra_mode" in h5f:
+            scan["zebra_mode"] = h5f["/entry1/zebra_mode"][0].decode()
+        else:
+            scan["zebra_mode"] = "nb"
+
+        # overwrite zebra_mode from cami
+        if cami_meta is not None:
+            if "zebra_mode" in cami_meta:
+                scan["zebra_mode"] = cami_meta["zebra_mode"][0]
+
+        if "/entry1/control/Monitor" in h5f:
+            scan["monitor"] = h5f["/entry1/control/Monitor"][0]
+        else:  # old path
+            scan["monitor"] = h5f["/entry1/control/data"][0]
+
+        scan["idx"] = 1
+
+        if "/entry1/sample/rotation_angle" in h5f:
+            scan["omega"] = h5f["/entry1/sample/rotation_angle"][:]
+        else:
+            scan["omega"] = h5f["/entry1/area_detector2/rotation_angle"][:]
+        if len(scan["omega"]) == 1:
+            scan["omega"] = np.ones(n) * scan["omega"]
+
+        scan["gamma"] = h5f["/entry1/ZEBRA/area_detector2/polar_angle"][:]
+        scan["twotheta"] = h5f["/entry1/ZEBRA/area_detector2/polar_angle"][:]
+        if len(scan["gamma"]) == 1:
+            scan["gamma"] = np.ones(n) * scan["gamma"]
+            scan["twotheta"] = np.ones(n) * scan["twotheta"]
+        scan["nu"] = h5f["/entry1/ZEBRA/area_detector2/tilt_angle"][:1]
+        scan["ddist"] = h5f["/entry1/ZEBRA/area_detector2/distance"][:1]
+        scan["wave"] = h5f["/entry1/ZEBRA/monochromator/wavelength"][:1]
+        scan["chi"] = h5f["/entry1/sample/chi"][:]
+        if len(scan["chi"]) == 1:
+            scan["chi"] = np.ones(n) * scan["chi"]
+        scan["phi"] = h5f["/entry1/sample/phi"][:]
+        if len(scan["phi"]) == 1:
+            scan["phi"] = np.ones(n) * scan["phi"]
+        scan["ub"] = h5f["/entry1/sample/UB"][:].reshape(3, 3)
+        scan["name"] = h5f["/entry1/sample/name"][0].decode()
+        scan["cell"] = h5f["/entry1/sample/cell"][:]
+
+        if n == 1:
+            # a default motor for a single frame file
+            scan["scan_motor"] = "omega"
+        else:
+            for var in ("omega", "gamma", "nu", "chi", "phi"):
+                if abs(scan[var][0] - scan[var][-1]) > 0.1:
+                    scan["scan_motor"] = var
+                    break
+            else:
+                raise ValueError("No angles that vary")
+
+        scan["scan_motors"] = [scan["scan_motor"], ]
 
         # optional parameters
         if "/entry1/sample/magnetic_field" in h5f:
-            det_data["magnetic_field"] = h5f["/entry1/sample/magnetic_field"][:]
+            scan["mf"] = h5f["/entry1/sample/magnetic_field"][:]
 
         if "/entry1/sample/temperature" in h5f:
-            det_data["temperature"] = h5f["/entry1/sample/temperature"][:]
+            scan["temp"] = h5f["/entry1/sample/temperature"][:]
+        elif "/entry1/sample/Ts/value" in h5f:
+            scan["temp"] = h5f["/entry1/sample/Ts/value"][:]
 
-    return det_data
+        # overwrite metadata from .cami
+        if cami_meta is not None:
+            if "crystal" in cami_meta:
+                cami_meta_crystal = cami_meta["crystal"]
+                if "name" in cami_meta_crystal:
+                    scan["name"] = cami_meta_crystal["name"]
+                if "UB" in cami_meta_crystal:
+                    scan["ub"] = cami_meta_crystal["UB"]
+                if "cell" in cami_meta_crystal:
+                    scan["cell"] = cami_meta_crystal["cell"]
+                if "lambda" in cami_meta_crystal:
+                    scan["wave"] = cami_meta_crystal["lambda"]
+
+            if "detector parameters" in cami_meta:
+                cami_meta_detparam = cami_meta["detector parameters"]
+                if "dist2" in cami_meta_detparam:
+                    scan["ddist"] = cami_meta_detparam["dist2"]
+
+    return scan
+
+
+def fit_event(scan, fr_from, fr_to, y_from, y_to, x_from, x_to):
+    data_roi = scan["counts"][fr_from:fr_to, y_from:y_to, x_from:x_to]
+
+    model = GaussianModel()
+    fr = np.arange(fr_from, fr_to)
+    counts_per_fr = np.sum(data_roi, axis=(1, 2))
+    params = model.guess(counts_per_fr, fr)
+    result = model.fit(counts_per_fr, x=fr, params=params)
+    frC = result.params["center"].value
+    intensity = result.params["height"].value
+
+    counts_std = counts_per_fr.std()
+    counts_mean = counts_per_fr.mean()
+    snr = 0 if counts_std == 0 else counts_mean / counts_std
+
+    model = Gaussian2dModel()
+    xs, ys = np.meshgrid(np.arange(x_from, x_to), np.arange(y_from, y_to))
+    xs = xs.flatten()
+    ys = ys.flatten()
+    counts = np.sum(data_roi, axis=0).flatten()
+    params = model.guess(counts, xs, ys)
+    result = model.fit(counts, x=xs, y=ys, params=params)
+    xC = result.params["centerx"].value
+    yC = result.params["centery"].value
+
+    scan["fit"] = {"frame": frC, "x_pos": xC, "y_pos": yC, "intensity": intensity, "snr": snr}

pyzebra/param_study_moduls.py

@@ -1,488 +0,0 @@
-import pickle
-
-import matplotlib as mpl
-import matplotlib.pyplot as plt
-import numpy as np
-import pandas as pd
-import scipy.io as sio
-import uncertainties as u
-from mpl_toolkits.mplot3d import Axes3D  # dont delete, otherwise waterfall wont work
-import collections
-
-from .ccl_io import load_1D
-
-def create_tuples(x, y, y_err):
-    """creates tuples for sorting and merginng of the data
-    Counts need to be normalized to monitor before"""
-    t = list()
-    for i in range(len(x)):
-        tup = (x[i], y[i], y_err[i])
-        t.append(tup)
-    return t
-
-
-def load_dats(filepath):
-    """reads the txt file, get headers and data
-    :arg filepath to txt file or list of filepaths to the files
-    :return ccl like dictionary"""
-    if isinstance(filepath, str):
-        data_type = "txt"
-        file_list = list()
-        with open(filepath, "r") as infile:
-            col_names = next(infile).split(",")
-            col_names = [col_names[i].rstrip() for i in range(len(col_names))]
-            for line in infile:
-                if "END" in line:
-                    break
-                file_list.append(tuple(line.split(",")))
-    elif isinstance(filepath, list):
-        data_type = "list"
-        file_list = filepath
-    dict1 = {}
-    for i in range(len(file_list)):
-        if not dict1:
-            if data_type == "txt":
-                dict1 = load_1D(file_list[0][0])
-            else:
-                dict1 = load_1D(file_list[0])
-        else:
-            if data_type == "txt":
-                dict1 = add_dict(dict1, load_1D(file_list[i][0]))
-            else:
-
-                dict1 = add_dict(dict1, load_1D(file_list[i]))
-        dict1["scan"][i + 1]["params"] = {}
-        if data_type == "txt":
-            for x in range(len(col_names) - 1):
-                dict1["scan"][i + 1]["params"][col_names[x + 1]] = float(file_list[i][x + 1])
-    return dict1
-
-
-def create_dataframe(dict1, variables):
-    """Creates pandas dataframe from the dictionary
-    :arg ccl like dictionary
-    :return pandas dataframe"""
-    # create dictionary to which we pull only wanted items before transforming it to pd.dataframe
-    pull_dict = {}
-    pull_dict["filenames"] = list()
-    for keys in variables:
-        for item in variables[keys]:
-            pull_dict[item] = list()
-    pull_dict["fit_area"] = list()
-    pull_dict["int_area"] = list()
-    pull_dict["Counts"] = list()
-
-    for keys in pull_dict:
-        print(keys)
-
-    # populate the dict
-    for keys in dict1["scan"]:
-        if "file_of_origin" in dict1["scan"][keys]:
-            pull_dict["filenames"].append(dict1["scan"][keys]["file_of_origin"].split("/")[-1])
-        else:
-            pull_dict["filenames"].append(dict1["meta"]["original_filename"].split("/")[-1])
-
-        pull_dict["fit_area"].append(dict1["scan"][keys]["fit"]["fit_area"])
-        pull_dict["int_area"].append(dict1["scan"][keys]["fit"]["int_area"])
-        pull_dict["Counts"].append(dict1["scan"][keys]["Counts"])
-        for key in variables:
-            for i in variables[key]:
-                pull_dict[i].append(_finditem(dict1["scan"][keys], i))
-
-    return pd.DataFrame(data=pull_dict)
-
-
-def sort_dataframe(dataframe, sorting_parameter):
-    """sorts the data frame and resets index"""
-    data = dataframe.sort_values(by=sorting_parameter)
-    data = data.reset_index(drop=True)
-    return data
-
-
-def make_graph(data, sorting_parameter, style):
-    """Makes the graph from the data based on style and sorting parameter
-    :arg data : pandas dataframe with data after sorting
-    :arg sorting_parameter to pull the correct variable and name
-    :arg style of the graph - waterfall, scatter, heatmap
-    :return matplotlib figure"""
-    if style == "waterfall":
-        mpl.rcParams["legend.fontsize"] = 10
-        fig = plt.figure()
-        ax = fig.gca(projection="3d")
-        for i in range(len(data)):
-            x = data["om"][i]
-            z = data["Counts"][i]
-            yy = [data[sorting_parameter][i]] * len(x)
-            ax.plot(x, yy, z, label=str("%s = %f" % (sorting_parameter, yy[i])))
-
-        ax.legend()
-        ax.set_xlabel("Omega")
-        ax.set_ylabel(sorting_parameter)
-        ax.set_zlabel("counts")
-
-    elif style == "scatter":
-        fig = plt.figure()
-        plt.errorbar(
-            data[sorting_parameter],
-            [data["fit_area"][i].n for i in range(len(data["fit_area"]))],
-            [data["fit_area"][i].s for i in range(len(data["fit_area"]))],
-            capsize=5,
-            ecolor="green",
-        )
-        plt.xlabel(str(sorting_parameter))
-        plt.ylabel("Intesity")
-
-    elif style == "heat":
-        new_om = list()
-        for i in range(len(data)):
-            new_om = np.append(new_om, np.around(data["om"][i], 2), axis=0)
-        unique_om = np.unique(new_om)
-        color_matrix = np.zeros(shape=(len(data), len(unique_om)))
-        for i in range(len(data)):
-            for j in range(len(data["om"][i])):
-                if np.around(data["om"][i][j], 2) in np.unique(new_om):
-                    color_matrix[i, j] = data["Counts"][i][j]
-                else:
-                    continue
-
-        fig = plt.figure()
-        plt.pcolormesh(unique_om, data[sorting_parameter], color_matrix, shading="gouraud")
-        plt.xlabel("omega")
-        plt.ylabel(sorting_parameter)
-        plt.colorbar()
-        plt.clim(color_matrix.mean(), color_matrix.max())
-
-    return fig
-
-
-def save_dict(obj, name):
-    """saves dictionary as pickle file in binary format
-    :arg obj - object to save
-    :arg name - name of the file
-    NOTE: path should be added later"""
-    with open(name + ".pkl", "wb") as f:
-        pickle.dump(obj, f, pickle.HIGHEST_PROTOCOL)
-
-
-def load_dict(name):
-    """load dictionary from picle file
-    :arg name - name of the file to load
-    NOTE: expect the file in the same folder, path should be added later
-    :return dictionary"""
-    with open(name + ".pkl", "rb") as f:
-        return pickle.load(f)
-
-
-# pickle, mat, h5, txt, csv, json
-def save_table(data, filetype, name, path=None):
-    print("Saving: ", filetype)
-    path = "" if path is None else path
-    if filetype == "pickle":
-        # to work with uncertanities, see uncertanity module
-        with open(path + name + ".pkl", "wb") as f:
-            pickle.dump(data, f, pickle.HIGHEST_PROTOCOL)
-    if filetype == "mat":
-        # matlab doesent allow some special character to be in var names, also cant start with
-        # numbers, in need, add some to the romove_character list
-        data["fit_area_nom"] = [data["fit_area"][i].n for i in range(len(data["fit_area"]))]
-        data["fit_area_err"] = [data["fit_area"][i].s for i in range(len(data["fit_area"]))]
-        data["int_area_nom"] = [data["int_area"][i].n for i in range(len(data["int_area"]))]
-        data["int_area_err"] = [data["int_area"][i].s for i in range(len(data["int_area"]))]
-        data = data.drop(columns=["fit_area", "int_area"])
-        remove_characters = [" ", "[", "]", "{", "}", "(", ")"]
-        for character in remove_characters:
-            data.columns = [
-                data.columns[i].replace(character, "") for i in range(len(data.columns))
-            ]
-        sio.savemat((path + name + ".mat"), {name: col.values for name, col in data.items()})
-    if filetype == "csv" or "txt":
-        data["fit_area_nom"] = [data["fit_area"][i].n for i in range(len(data["fit_area"]))]
-        data["fit_area_err"] = [data["fit_area"][i].s for i in range(len(data["fit_area"]))]
-        data["int_area_nom"] = [data["int_area"][i].n for i in range(len(data["int_area"]))]
-        data["int_area_err"] = [data["int_area"][i].s for i in range(len(data["int_area"]))]
-        data = data.drop(columns=["fit_area", "int_area", "om", "Counts"])
-        if filetype == "csv":
-            data.to_csv(path + name + ".csv")
-        if filetype == "txt":
-            with open((path + name + ".txt"), "w") as outfile:
-                data.to_string(outfile)
-    if filetype == "h5":
-        hdf = pd.HDFStore((path + name + ".h5"))
-        hdf.put("data", data)
-        hdf.close()
-    if filetype == "json":
-        data.to_json((path + name + ".json"))
-
-
-def normalize(scan, monitor):
-    """Normalizes the measurement to monitor, checks if sigma exists, otherwise creates it
-    :arg dict : dictionary to from which to tkae the scan
-    :arg key : which scan to normalize from dict1
-    :arg monitor : final monitor
-    :return counts - normalized counts
-    :return sigma - normalized sigma"""
-
-    counts = np.array(scan["Counts"])
-    sigma = np.sqrt(counts) if "sigma" not in scan else scan["sigma"]
-    monitor_ratio = monitor / scan["monitor"]
-    scaled_counts = counts * monitor_ratio
-    scaled_sigma = np.array(sigma) * monitor_ratio
-
-    return scaled_counts, scaled_sigma
-
-
-def merge(scan1, scan2, keep=True, monitor=100000):
-    """merges the two tuples and sorts them, if om value is same, Counts value is average
-    averaging is propagated into sigma if dict1 == dict2, key[1] is deleted after merging
-    :arg dict1 : dictionary to which measurement will be merged
-    :arg dict2 : dictionary from which measurement will be merged
-    :arg scand_dict_result : result of scan_dict after auto function
-    :arg keep : if true, when monitors are same, does not change it, if flase, takes monitor
-    always
-    :arg monitor : final monitor after merging
-    note: dict1 and dict2 can be same dict
-    :return dict1 with merged scan"""
-
-    if keep:
-        if scan1["monitor"] == scan2["monitor"]:
-            monitor = scan1["monitor"]
-
-    # load om and Counts
-    x1, x2 = scan1["om"], scan2["om"]
-    cor_y1, y_err1 = normalize(scan1, monitor=monitor)
-    cor_y2, y_err2 = normalize(scan2, monitor=monitor)
-    # creates touples (om, Counts, sigma) for sorting and further processing
-    tuple_list = create_tuples(x1, cor_y1, y_err1) + create_tuples(x2, cor_y2, y_err2)
-    # Sort the list on om and add 0 0 0 tuple to the last position
-    sorted_t = sorted(tuple_list, key=lambda tup: tup[0])
-    sorted_t.append((0, 0, 0))
-    om, Counts, sigma = [], [], []
-    seen = list()
-    for i in range(len(sorted_t) - 1):
-        if sorted_t[i][0] not in seen:
-            if sorted_t[i][0] != sorted_t[i + 1][0]:
-                om = np.append(om, sorted_t[i][0])
-                Counts = np.append(Counts, sorted_t[i][1])
-                sigma = np.append(sigma, sorted_t[i][2])
-            else:
-                om = np.append(om, sorted_t[i][0])
-                counts1, counts2 = sorted_t[i][1], sorted_t[i + 1][1]
-                sigma1, sigma2 = sorted_t[i][2], sorted_t[i + 1][2]
-                count_err1 = u.ufloat(counts1, sigma1)
-                count_err2 = u.ufloat(counts2, sigma2)
-                avg = (count_err1 + count_err2) / 2
-                Counts = np.append(Counts, avg.n)
-                sigma = np.append(sigma, avg.s)
-                seen.append(sorted_t[i][0])
-        else:
-            continue
-    scan1["om"] = om
-    scan1["Counts"] = Counts
-    scan1["sigma"] = sigma
-    scan1["monitor"] = monitor
-    print("merging done")
-
-
-def add_dict(dict1, dict2):
-    """adds two dictionaries, meta of the new is saved as meata+original_filename and
-    measurements are shifted to continue with numbering of first dict
-    :arg dict1 : dictionarry to add to
-    :arg dict2 : dictionarry from which to take the measurements
-    :return dict1 : combined dictionary
-    Note: dict1 must be made from ccl, otherwise we would have to change the structure of loaded
-    dat file"""
-    try:
-        if dict1["meta"]["zebra_mode"] != dict2["meta"]["zebra_mode"]:
-            print("You are trying to add scans measured with different zebra modes")
-            return
-    # this is for the qscan case
-    except KeyError:
-        print("Zebra mode not specified")
-    max_measurement_dict1 = max([keys for keys in dict1["scan"]])
-    new_filenames = np.arange(
-        max_measurement_dict1 + 1, max_measurement_dict1 + 1 + len(dict2["scan"])
-    )
-    new_meta_name = "meta" + str(dict2["meta"]["original_filename"])
-    if new_meta_name not in dict1:
-        for keys, name in zip(dict2["scan"], new_filenames):
-            dict2["scan"][keys]["file_of_origin"] = str(dict2["meta"]["original_filename"])
-            dict1["scan"][name] = dict2["scan"][keys]
-
-        dict1[new_meta_name] = dict2["meta"]
-    else:
-        raise KeyError(
-            str(
-                "The file %s has alredy been added to %s"
-                % (dict2["meta"]["original_filename"], dict1["meta"]["original_filename"])
-            )
-        )
-    return dict1
-
-
-def auto(dict):
-    """takes just unique tuples from all tuples in dictionary returend by scan_dict
-    intendet for automatic merge if you doesent want to specify what scans to merge together
-    args: dict - dictionary from scan_dict function
-    :return dict - dict without repetitions"""
-    for keys in dict:
-        tuple_list = dict[keys]
-        new = list()
-        for i in range(len(tuple_list)):
-            if tuple_list[0][0] == tuple_list[i][0]:
-                new.append(tuple_list[i])
-        dict[keys] = new
-    return dict
-
-
-def scan_dict(dict, precision=0.5):
-    """scans dictionary for duplicate angles indexes
-    :arg dict : dictionary to scan
-    :arg precision : in deg, sometimes angles are zero so its easier this way, instead of
-    checking zero division
-    :return  dictionary with matching scans, if there are none, the dict is empty
-    note: can be checked by "not d", true if empty
-    """
-
-    if dict["meta"]["zebra_mode"] == "bi":
-        angles = ["twotheta_angle", "omega_angle", "chi_angle", "phi_angle"]
-    elif dict["meta"]["zebra_mode"] == "nb":
-        angles = ["gamma_angle", "omega_angle", "nu_angle"]
-    else:
-        print("Unknown zebra mode")
-        return
-
-    d = {}
-    for i in dict["scan"]:
-        for j in dict["scan"]:
-            if dict["scan"][i] != dict["scan"][j]:
-                itup = list()
-                for k in angles:
-                    itup.append(abs(abs(dict["scan"][i][k]) - abs(dict["scan"][j][k])))
-
-                if all(i <= precision for i in itup):
-                    print(itup)
-                    print([dict["scan"][i][k] for k in angles])
-                    print([dict["scan"][j][k] for k in angles])
-                    if str([np.around(dict["scan"][i][k], 0) for k in angles]) not in d:
-                        d[str([np.around(dict["scan"][i][k], 0) for k in angles])] = list()
-                        d[str([np.around(dict["scan"][i][k], 0) for k in angles])].append((i, j))
-                    else:
-                        d[str([np.around(dict["scan"][i][k], 0) for k in angles])].append((i, j))
-
-                else:
-                    pass
-
-            else:
-                continue
-
-    return d
-
-
-def _finditem(obj, key):
-    if key in obj:
-        return obj[key]
-    for k, v in obj.items():
-        if isinstance(v, dict):
-            item = _finditem(v, key)
-            if item is not None:
-                return item
-
-
-def most_common(lst):
-    return max(set(lst), key=lst.count)
-
-
-def variables(dictionary):
-    """Funcrion to guess what variables will be used in the param study
-    i call pripary variable the one the array like variable, usually omega
-    and secondary the slicing variable, different for each scan,for example temperature"""
-    # find all variables that are in all scans
-    stdev_precision = 0.05
-    all_vars = list()
-    for keys in dictionary["scan"]:
-        all_vars.append([key for key in dictionary["scan"][keys] if key != "params"])
-        if dictionary["scan"][keys]["params"]:
-            all_vars.append(key for key in dictionary["scan"][keys]["params"])
-
-    all_vars = [i for sublist in all_vars for i in sublist]
-    # get the ones that are in all scans
-    b = collections.Counter(all_vars)
-    inall = [key for key in b if b[key] == len(dictionary["scan"])]
-    # delete those that are obviously wrong
-    wrong = [
-        "NP",
-        "Counts",
-        "Monitor1",
-        "Monitor2",
-        "Monitor3",
-        "h_index",
-        "l_index",
-        "k_index",
-        "n_points",
-        "monitor",
-        "Time",
-        "omega_angle",
-        "twotheta_angle",
-        "chi_angle",
-        "phi_angle",
-        "nu_angle",
-    ]
-    inall_red = [i for i in inall if i not in wrong]
-
-    # check for primary variable, needs to be list, we dont suspect the
-    # primary variable be as a parameter (be in scan[params])
-    primary_candidates = list()
-    for key in dictionary["scan"]:
-        for i in inall_red:
-            if isinstance(_finditem(dictionary["scan"][key], i), list):
-                if np.std(_finditem(dictionary["scan"][key], i)) > stdev_precision:
-                    primary_candidates.append(i)
-    # check which of the primary are in every scan
-    primary_candidates = collections.Counter(primary_candidates)
-    second_round_primary_candidates = [
-        key for key in primary_candidates if primary_candidates[key] == len(dictionary["scan"])
-    ]
-
-    if len(second_round_primary_candidates) == 1:
-        print("We've got a primary winner!", second_round_primary_candidates)
-    else:
-        print("Still not sure with primary:(", second_round_primary_candidates)
-
-    # check for secondary variable, we suspect a float\int or not changing array
-    # we dont need to check for primary ones
-    secondary_candidates = [i for i in inall_red if i not in second_round_primary_candidates]
-    # print("secondary candidates", secondary_candidates)
-    # select arrays and floats and ints
-    second_round_secondary_candidates = list()
-    for key in dictionary["scan"]:
-        for i in secondary_candidates:
-            if isinstance(_finditem(dictionary["scan"][key], i), float):
-                second_round_secondary_candidates.append(i)
-            elif isinstance(_finditem(dictionary["scan"][key], i), int):
-                second_round_secondary_candidates.append(i)
-            elif isinstance(_finditem(dictionary["scan"][key], i), list):
-                if np.std(_finditem(dictionary["scan"][key], i)) < stdev_precision:
-                    second_round_secondary_candidates.append(i)
-
-    second_round_secondary_candidates = collections.Counter(second_round_secondary_candidates)
-    second_round_secondary_candidates = [
-        key
-        for key in second_round_secondary_candidates
-        if second_round_secondary_candidates[key] == len(dictionary["scan"])
-    ]
-    # print("secondary candidates after second round", second_round_secondary_candidates)
-    # now we check if they vary between the scans
-    third_round_sec_candidates = list()
-    for i in second_round_secondary_candidates:
-        check_array = list()
-        for keys in dictionary["scan"]:
-            check_array.append(np.average(_finditem(dictionary["scan"][keys], i)))
-        # print(i, check_array, np.std(check_array))
-        if np.std(check_array) > stdev_precision:
-            third_round_sec_candidates.append(i)
-    if len(third_round_sec_candidates) == 1:
-        print("We've got a secondary winner!", third_round_sec_candidates)
-    else:
-        print("Still not sure with secondary :(", third_round_sec_candidates)
-
-    return {"primary": second_round_primary_candidates, "secondary": third_round_sec_candidates}

pyzebra/sxtal_refgen.py

@@ -0,0 +1,483 @@
+import io
+import os
+import subprocess
+import tempfile
+from math import ceil, floor
+
+import numpy as np
+
+SXTAL_REFGEN_PATH = "/afs/psi.ch/project/sinq/rhel7/bin/Sxtal_Refgen"
+
+_zebraBI_default_geom = """GEOM          2 Bissecting - HiCHI
+BLFR          z-up
+DIST_UNITS    mm
+ANGL_UNITS    deg
+DET_TYPE      Point  ipsd 1
+DIST_DET      488
+DIM_XY         1.0   1.0    1    1
+GAPS_DET       0  0
+
+SETTING       1 0 0   0 1 0   0 0 1
+NUM_ANG  4
+ANG_LIMITS         Min      Max    Offset
+        Gamma     0.0    128.0     0.00
+        Omega     0.0     64.0     0.00
+        Chi      80.0    211.0     0.00
+        Phi       0.0    360.0     0.00
+
+DET_OFF      0       0      0
+"""
+
+_zebraNB_default_geom = """GEOM          3 Normal Beam
+BLFR          z-up
+DIST_UNITS    mm
+ANGL_UNITS    deg
+DET_TYPE      Point  ipsd 1
+DIST_DET      448
+DIM_XY         1.0   1.0    1    1
+GAPS_DET       0  0
+
+SETTING       1 0 0   0 1 0   0 0 1
+NUM_ANG  3
+ANG_LIMITS         Min      Max    Offset
+        Gamma      0.0     128.0     0.00
+        Omega   -180.0     180.0     0.00
+        Nu       -15.0      15.0     0.00
+
+DET_OFF      0       0      0
+"""
+
+_zebra_default_cfl = """TITLE mymaterial
+SPGR  P 63 2 2
+CELL  5.73 5.73 11.89 90 90 120
+
+WAVE 1.383
+
+UBMAT
+0.000000    0.000000    0.084104
+0.000000    0.174520   -0.000000
+0.201518    0.100759    0.000000
+
+INSTR  zebra.geom
+
+ORDER   1 2 3
+
+ANGOR   gamma
+
+HLIM    -25 25 -25 25 -25 25
+SRANG   0.0  0.7
+
+Mag_Structure
+lattiCE P 1
+kvect        0.0 0.0 0.0
+magcent
+symm  x,y,z
+msym  u,v,w, 0.0
+End_Mag_Structure
+"""
+
+
+def get_zebraBI_default_geom_file():
+    return io.StringIO(_zebraBI_default_geom)
+
+
+def get_zebraNB_default_geom_file():
+    return io.StringIO(_zebraNB_default_geom)
+
+
+def get_zebra_default_cfl_file():
+    return io.StringIO(_zebra_default_cfl)
+
+
+def read_geom_file(fileobj):
+    ang_lims = dict()
+    for line in fileobj:
+        if "!" in line:  # remove comments that start with ! sign
+            line, _ = line.split(sep="!", maxsplit=1)
+
+        if line.startswith("GEOM"):
+            _, val = line.split(maxsplit=1)
+            if val.startswith("2"):
+                ang_lims["geom"] = "bi"
+            else:  # val.startswith("3")
+                ang_lims["geom"] = "nb"
+
+        elif line.startswith("ANG_LIMITS"):
+            # read angular limits
+            for line in fileobj:
+                if not line or line.isspace():
+                    break
+
+                ang, ang_min, ang_max, ang_offset = line.split()
+                ang_lims[ang.lower()] = [ang_min, ang_max, ang_offset]
+
+            if "2theta" in ang_lims:  # treat 2theta as gamma
+                ang_lims["gamma"] = ang_lims.pop("2theta")
+
+    return ang_lims
+
+
+def export_geom_file(path, ang_lims, template=None):
+    if ang_lims["geom"] == "bi":
+        template_file = get_zebraBI_default_geom_file()
+        n_ang = 4
+    else:  # ang_lims["geom"] == "nb"
+        template_file = get_zebraNB_default_geom_file()
+        n_ang = 3
+
+    if template is not None:
+        template_file = template
+
+    with open(path, "w") as out_file:
+        for line in template_file:
+            out_file.write(line)
+
+            if line.startswith("ANG_LIMITS"):
+                for _ in range(n_ang):
+                    next_line = next(template_file)
+                    ang, _, _, _ = next_line.split()
+
+                    if ang == "2theta":  # treat 2theta as gamma
+                        ang = "Gamma"
+                    vals = ang_lims[ang.lower()]
+
+                    out_file.write(f"{'':<8}{ang:<10}{vals[0]:<10}{vals[1]:<10}{vals[2]:<10}\n")
+
+
+def calc_ub_matrix(params):
+    with tempfile.TemporaryDirectory() as temp_dir:
+        cfl_file = os.path.join(temp_dir, "ub_matrix.cfl")
+
+        with open(cfl_file, "w") as fileobj:
+            for key, value in params.items():
+                fileobj.write(f"{key} {value}\n")
+
+        comp_proc = subprocess.run(
+            [SXTAL_REFGEN_PATH, cfl_file],
+            cwd=temp_dir,
+            check=True,
+            stdout=subprocess.PIPE,
+            stderr=subprocess.STDOUT,
+            text=True,
+        )
+        print(" ".join(comp_proc.args))
+        print(comp_proc.stdout)
+
+        sfa_file = os.path.join(temp_dir, "ub_matrix.sfa")
+        ub_matrix = []
+        with open(sfa_file, "r") as fileobj:
+            for line in fileobj:
+                if "BL_M" in line:  # next 3 lines contain the matrix
+                    for _ in range(3):
+                        next_line = next(fileobj)
+                        *vals, _ = next_line.split(maxsplit=3)
+                        ub_matrix.extend(vals)
+
+    return ub_matrix
+
+
+def read_cfl_file(fileobj):
+    params = {
+        "SPGR": None,
+        "CELL": None,
+        "WAVE": None,
+        "UBMAT": None,
+        "HLIM": None,
+        "SRANG": None,
+        "lattiCE": None,
+        "kvect": None,
+    }
+    param_names = tuple(params)
+
+    for line in fileobj:
+        line = line.strip()
+        if "!" in line:  # remove comments that start with ! sign
+            line, _ = line.split(sep="!", maxsplit=1)
+
+        if line.startswith(param_names):
+            if line.startswith("UBMAT"):  # next 3 lines contain the matrix
+                param, val = "UBMAT", []
+                for _ in range(3):
+                    next_line = next(fileobj).strip()
+                    val.extend(next_line.split(maxsplit=2))
+            else:
+                param, val = line.split(maxsplit=1)
+
+            params[param] = val
+
+    return params
+
+
+def read_cif_file(fileobj):
+    params = {"SPGR": None, "CELL": None, "ATOM": []}
+
+    cell_params = {
+        "_cell_length_a": None,
+        "_cell_length_b": None,
+        "_cell_length_c": None,
+        "_cell_angle_alpha": None,
+        "_cell_angle_beta": None,
+        "_cell_angle_gamma": None,
+    }
+    cell_param_names = tuple(cell_params)
+
+    atom_param_pos = {
+        "_atom_site_label": 0,
+        "_atom_site_type_symbol": None,
+        "_atom_site_fract_x": None,
+        "_atom_site_fract_y": None,
+        "_atom_site_fract_z": None,
+        "_atom_site_U_iso_or_equiv": None,
+        "_atom_site_occupancy": None,
+    }
+    atom_param_names = tuple(atom_param_pos)
+
+    for line in fileobj:
+        line = line.strip()
+        if line.startswith("_space_group_name_H-M_alt"):
+            _, val = line.split(maxsplit=1)
+            params["SPGR"] = val.strip("'")
+
+        elif line.startswith(cell_param_names):
+            param, val = line.split(maxsplit=1)
+            cell_params[param] = val
+
+        elif line.startswith("_atom_site_label"):  # assume this is the start of atom data
+            for ind, line in enumerate(fileobj, start=1):
+                line = line.strip()
+
+                # read fields
+                if line.startswith("_atom_site"):
+                    if line.startswith(atom_param_names):
+                        atom_param_pos[line] = ind
+                    continue
+
+                # read data till an empty line
+                if not line:
+                    break
+                vals = line.split()
+                params["ATOM"].append(" ".join([vals[ind] for ind in atom_param_pos.values()]))
+
+    if None not in cell_params.values():
+        params["CELL"] = " ".join(cell_params.values())
+
+    return params
+
+
+def export_cfl_file(path, params, template=None):
+    param_names = tuple(params)
+    if template is None:
+        template_file = get_zebra_default_cfl_file()
+    else:
+        template_file = template
+
+    atom_done = False
+    with open(path, "w") as out_file:
+        for line in template_file:
+            if line.startswith(param_names):
+                if line.startswith("UBMAT"):  # only UBMAT values are not on the same line
+                    out_file.write(line)
+                    for i in range(3):
+                        next(template_file)
+                        out_file.write(" ".join(params["UBMAT"][3 * i : 3 * (i + 1)]) + "\n")
+
+                elif line.startswith("ATOM"):
+                    if "ATOM" in params:
+                        # replace all ATOM with values in params
+                        while line.startswith("ATOM"):
+                            line = next(template_file)
+                        for atom_line in params["ATOM"]:
+                            out_file.write(f"ATOM {atom_line}\n")
+                        atom_done = True
+
+                else:
+                    param, _ = line.split(maxsplit=1)
+                    out_file.write(f"{param} {params[param]}\n")
+
+            elif line.startswith("INSTR"):
+                # replace it with a default name
+                out_file.write("INSTR  zebra.geom\n")
+
+            else:
+                out_file.write(line)
+
+        # append ATOM data if it's present and a template did not contain it
+        if "ATOM" in params and not atom_done:
+            out_file.write("\n")
+            for atom_line in params["ATOM"]:
+                out_file.write(f"ATOM {atom_line}\n")
+
+
+def sort_hkl_file_bi(file_in, file_out, priority, chunks):
+    with open(file_in) as fileobj:
+        file_in_data = fileobj.readlines()
+
+    data = np.genfromtxt(file_in, skip_header=3)
+    stt = data[:, 4]
+    omega = data[:, 5]
+    chi = data[:, 6]
+    phi = data[:, 7]
+
+    lines = file_in_data[3:]
+    lines_update = []
+
+    angles = {"2theta": stt, "omega": omega, "chi": chi, "phi": phi}
+
+    # Reverse flag
+    to_reverse = False
+    to_reverse_p2 = False
+    to_reverse_p3 = False
+
+    # Get indices within first priority
+    ang_p1 = angles[priority[0]]
+    begin_p1 = floor(min(ang_p1))
+    end_p1 = ceil(max(ang_p1))
+    delta_p1 = chunks[0]
+    for p1 in range(begin_p1, end_p1, delta_p1):
+        ind_p1 = [j for j, x in enumerate(ang_p1) if p1 <= x and x < p1 + delta_p1]
+
+        stt_new = [stt[x] for x in ind_p1]
+        omega_new = [omega[x] for x in ind_p1]
+        chi_new = [chi[x] for x in ind_p1]
+        phi_new = [phi[x] for x in ind_p1]
+        lines_new = [lines[x] for x in ind_p1]
+
+        angles_p2 = {"stt": stt_new, "omega": omega_new, "chi": chi_new, "phi": phi_new}
+
+        # Get indices for second priority
+        ang_p2 = angles_p2[priority[1]]
+        if len(ang_p2) > 0 and to_reverse_p2:
+            begin_p2 = ceil(max(ang_p2))
+            end_p2 = floor(min(ang_p2))
+            delta_p2 = -chunks[1]
+        elif len(ang_p2) > 0 and not to_reverse_p2:
+            end_p2 = ceil(max(ang_p2))
+            begin_p2 = floor(min(ang_p2))
+            delta_p2 = chunks[1]
+        else:
+            end_p2 = 0
+            begin_p2 = 0
+            delta_p2 = 1
+
+        to_reverse_p2 = not to_reverse_p2
+
+        for p2 in range(begin_p2, end_p2, delta_p2):
+            min_p2 = min([p2, p2 + delta_p2])
+            max_p2 = max([p2, p2 + delta_p2])
+            ind_p2 = [j for j, x in enumerate(ang_p2) if min_p2 <= x and x < max_p2]
+
+            stt_new2 = [stt_new[x] for x in ind_p2]
+            omega_new2 = [omega_new[x] for x in ind_p2]
+            chi_new2 = [chi_new[x] for x in ind_p2]
+            phi_new2 = [phi_new[x] for x in ind_p2]
+            lines_new2 = [lines_new[x] for x in ind_p2]
+
+            angles_p3 = {"stt": stt_new2, "omega": omega_new2, "chi": chi_new2, "phi": phi_new2}
+
+            # Get indices for third priority
+            ang_p3 = angles_p3[priority[2]]
+            if len(ang_p3) > 0 and to_reverse_p3:
+                begin_p3 = ceil(max(ang_p3)) + chunks[2]
+                end_p3 = floor(min(ang_p3)) - chunks[2]
+                delta_p3 = -chunks[2]
+            elif len(ang_p3) > 0 and not to_reverse_p3:
+                end_p3 = ceil(max(ang_p3)) + chunks[2]
+                begin_p3 = floor(min(ang_p3)) - chunks[2]
+                delta_p3 = chunks[2]
+            else:
+                end_p3 = 0
+                begin_p3 = 0
+                delta_p3 = 1
+
+            to_reverse_p3 = not to_reverse_p3
+
+            for p3 in range(begin_p3, end_p3, delta_p3):
+                min_p3 = min([p3, p3 + delta_p3])
+                max_p3 = max([p3, p3 + delta_p3])
+                ind_p3 = [j for j, x in enumerate(ang_p3) if min_p3 <= x and x < max_p3]
+
+                angle_new3 = [angles_p3[priority[3]][x] for x in ind_p3]
+
+                ind_final = [x for _, x in sorted(zip(angle_new3, ind_p3), reverse=to_reverse)]
+
+                to_reverse = not to_reverse
+
+                for i in ind_final:
+                    lines_update.append(lines_new2[i])
+
+    with open(file_out, "w") as fileobj:
+        for _ in range(3):
+            fileobj.write(file_in_data.pop(0))
+
+        fileobj.writelines(lines_update)
+
+
+def sort_hkl_file_nb(file_in, file_out, priority, chunks):
+    with open(file_in) as fileobj:
+        file_in_data = fileobj.readlines()
+
+    data = np.genfromtxt(file_in, skip_header=3)
+    gamma = data[:, 4]
+    omega = data[:, 5]
+    nu = data[:, 6]
+
+    lines = file_in_data[3:]
+    lines_update = []
+
+    angles = {"gamma": gamma, "omega": omega, "nu": nu}
+
+    to_reverse = False
+    to_reverse_p2 = False
+
+    # Get indices within first priority
+    ang_p1 = angles[priority[0]]
+    begin_p1 = floor(min(ang_p1))
+    end_p1 = ceil(max(ang_p1))
+    delta_p1 = chunks[0]
+    for p1 in range(begin_p1, end_p1, delta_p1):
+        ind_p1 = [j for j, x in enumerate(ang_p1) if p1 <= x and x < p1 + delta_p1]
+
+        # Get angles from within nu range
+        lines_new = [lines[x] for x in ind_p1]
+        gamma_new = [gamma[x] for x in ind_p1]
+        omega_new = [omega[x] for x in ind_p1]
+        nu_new = [nu[x] for x in ind_p1]
+
+        angles_p2 = {"gamma": gamma_new, "omega": omega_new, "nu": nu_new}
+
+        # Get indices for second priority
+        ang_p2 = angles_p2[priority[1]]
+        if len(gamma_new) > 0 and to_reverse_p2:
+            begin_p2 = ceil(max(ang_p2))
+            end_p2 = floor(min(ang_p2))
+            delta_p2 = -chunks[1]
+        elif len(gamma_new) > 0 and not to_reverse_p2:
+            end_p2 = ceil(max(ang_p2))
+            begin_p2 = floor(min(ang_p2))
+            delta_p2 = chunks[1]
+        else:
+            end_p2 = 0
+            begin_p2 = 0
+            delta_p2 = 1
+
+        to_reverse_p2 = not to_reverse_p2
+
+        for p2 in range(begin_p2, end_p2, delta_p2):
+            min_p2 = min([p2, p2 + delta_p2])
+            max_p2 = max([p2, p2 + delta_p2])
+            ind_p2 = [j for j, x in enumerate(ang_p2) if min_p2 <= x and x < max_p2]
+
+            angle_new2 = [angles_p2[priority[2]][x] for x in ind_p2]
+
+            ind_final = [x for _, x in sorted(zip(angle_new2, ind_p2), reverse=to_reverse)]
+
+            to_reverse = not to_reverse
+
+            for i in ind_final:
+                lines_update.append(lines_new[i])
+
+    with open(file_out, "w") as fileobj:
+        for _ in range(3):
+            fileobj.write(file_in_data.pop(0))
+
+        fileobj.writelines(lines_update)

pyzebra/utils.py

@@ -0,0 +1,17 @@
+import os
+
+SINQ_PATH = "/afs/psi.ch/project/sinqdata"
+ZEBRA_PROPOSALS_PATH = os.path.join(SINQ_PATH, "{year}/zebra/{proposal}")
+
+
+def find_proposal_path(proposal):
+    for entry in os.scandir(SINQ_PATH):
+        if entry.is_dir() and len(entry.name) == 4 and entry.name.isdigit():
+            proposal_path = ZEBRA_PROPOSALS_PATH.format(year=entry.name, proposal=proposal)
+            if os.path.isdir(proposal_path):
+                # found it
+                break
+    else:
+        raise ValueError(f"Can not find data for proposal '{proposal}'.")
+
+    return proposal_path

pyzebra/xtal.py

@@ -1,15 +1,5 @@
-import math
-
 import numpy as np
 from numba import njit
-from scipy.optimize import curve_fit
-
-import pyzebra
-
-try:
-    from matplotlib import pyplot as plt
-except ImportError:
-    print("matplotlib is not available")
 
 pi_r = 180 / np.pi
 
@@ -382,6 +372,27 @@ def ang2hkl(wave, ddist, gammad, om, ch, ph, nud, ub, x, y):
     return hkl
 
 
+def ang2hkl_1d(wave, ddist, ga, om, ch, ph, nu, ub):
+    """Calculate hkl-indices of a reflection from its position (angles) at the 1d-detector
+    """
+    z1 = z1frmd(wave, ga, om, ch, ph, nu)
+    ubinv = np.linalg.inv(ub)
+    hkl = ubinv @ z1
+
+    return hkl
+
+
+def ang_proc(wave, ddist, gammad, om, ch, ph, nud, x, y):
+    """Utility function to calculate ch, ph, ga, om
+    """
+    ga, nu = det2pol(ddist, gammad, nud, x, y)
+    z1 = z1frmd(wave, ga, om, ch, ph, nu)
+    ch2, ph2 = eqchph(z1)
+    ch, ph, ga, om = fixdnu(wave, z1, ch2, ph2, nu)
+
+    return ch, ph, ga, om
+
+
 def gauss(x, *p):
     """Defines Gaussian function
 
@@ -393,84 +404,3 @@ def gauss(x, *p):
     """
     A, mu, sigma = p
     return A * np.exp(-((x - mu) ** 2) / (2.0 * sigma ** 2))
-
-
-def box_int(file, box):
-    """Calculates center of the peak in the NB-geometry angles and Intensity of the peak
-
-    Args:
-        file name, box size [x0:xN, y0:yN, fr0:frN]
-
-    Returns:
-        gamma, omPeak, nu polar angles, Int and data for 3 fit plots
-    """
-
-    dat = pyzebra.read_detector_data(file)
-
-    sttC = dat["pol_angle"][0]
-    om = dat["rot_angle"]
-    nuC = dat["tlt_angle"][0]
-    ddist = dat["ddist"]
-
-    # defining indices
-    x0, xN, y0, yN, fr0, frN = box
-
-    # omega fit
-    om = dat["rot_angle"][fr0:frN]
-    cnts = np.sum(dat["data"][fr0:frN, y0:yN, x0:xN], axis=(1, 2))
-
-    p0 = [1.0, 0.0, 1.0]
-    coeff, var_matrix = curve_fit(gauss, range(len(cnts)), cnts, p0=p0)
-
-    frC = fr0 + coeff[1]
-    omF = dat["rot_angle"][math.floor(frC)]
-    omC = dat["rot_angle"][math.ceil(frC)]
-    frStep = frC - math.floor(frC)
-    omStep = omC - omF
-    omP = omF + omStep * frStep
-    Int = coeff[1] * abs(coeff[2] * omStep) * math.sqrt(2) * math.sqrt(np.pi)
-    # omega plot
-    x_fit = np.linspace(0, len(cnts), 100)
-    y_fit = gauss(x_fit, *coeff)
-    plt.figure()
-    plt.subplot(131)
-    plt.plot(range(len(cnts)), cnts)
-    plt.plot(x_fit, y_fit)
-    plt.ylabel("Intensity in the box")
-    plt.xlabel("Frame N of the box")
-    label = "om"
-    # gamma fit
-    sliceXY = dat["data"][fr0:frN, y0:yN, x0:xN]
-    sliceXZ = np.sum(sliceXY, axis=1)
-    sliceYZ = np.sum(sliceXY, axis=2)
-
-    projX = np.sum(sliceXZ, axis=0)
-    p0 = [1.0, 0.0, 1.0]
-    coeff, var_matrix = curve_fit(gauss, range(len(projX)), projX, p0=p0)
-    x = x0 + coeff[1]
-    # gamma plot
-    x_fit = np.linspace(0, len(projX), 100)
-    y_fit = gauss(x_fit, *coeff)
-    plt.subplot(132)
-    plt.plot(range(len(projX)), projX)
-    plt.plot(x_fit, y_fit)
-    plt.ylabel("Intensity in the box")
-    plt.xlabel("X-pixel of the box")
-
-    # nu fit
-    projY = np.sum(sliceYZ, axis=0)
-    p0 = [1.0, 0.0, 1.0]
-    coeff, var_matrix = curve_fit(gauss, range(len(projY)), projY, p0=p0)
-    y = y0 + coeff[1]
-    # nu plot
-    x_fit = np.linspace(0, len(projY), 100)
-    y_fit = gauss(x_fit, *coeff)
-    plt.subplot(133)
-    plt.plot(range(len(projY)), projY)
-    plt.plot(x_fit, y_fit)
-    plt.ylabel("Intensity in the box")
-    plt.xlabel("Y-pixel of the box")
-
-    ga, nu = pyzebra.det2pol(ddist, sttC, nuC, x, y)
-
-    return ga[0], omP, nu[0], Int

scripts/pyzebra-start.sh

@@ -0,0 +1,4 @@
+source /home/pyzebra/miniconda3/etc/profile.d/conda.sh
+
+conda activate prod
+pyzebra --port=80 --allow-websocket-origin=pyzebra.psi.ch:80 --spind-path=/home/pyzebra/spind

scripts/pyzebra-test-start.sh

@@ -0,0 +1,4 @@
+source /home/pyzebra/miniconda3/etc/profile.d/conda.sh
+
+conda activate test
+python ~/pyzebra/pyzebra/app/cli.py --allow-websocket-origin=pyzebra.psi.ch:5006 --spind-path=/home/pyzebra/spind

scripts/pyzebra-test.service

@@ -0,0 +1,11 @@
+[Unit]
+Description=pyzebra-test web server (runs on port 5006)
+
+[Service]
+Type=simple
+User=pyzebra
+ExecStart=/bin/bash /usr/local/sbin/pyzebra-test-start.sh
+Restart=always
+
+[Install]
+WantedBy=multi-user.target

scripts/pyzebra.service

@@ -0,0 +1,10 @@
+[Unit]
+Description=pyzebra web server
+
+[Service]
+Type=simple
+ExecStart=/bin/bash /usr/local/sbin/pyzebra-start.sh
+Restart=always
+
+[Install]
+WantedBy=multi-user.target