Error calculation action and bbq (#381)

* Tune error based on deviation of filtered BBQ data to the moving average
* Action error calculated from error on the spectral line

CHANGELOG.md

@@ -1,5 +1,13 @@
 # OMC3 Changelog
 
+#### 2022-11-08 - v0.7.0 - _jdilly_
+
+- Added:
+  - Tune error based on deviation of filtered BBQ data to the moving average
+    (over moving average window)
+  - Action error calculated from error on the spectral line
+    (which in turn is the same as NOISE)
+
 #### 2022-11-01 - v0.6.6
 
 - Bugfixes 

omc3/__init__.py

@@ -11,7 +11,7 @@
 __title__ = "omc3"
 __description__ = "An accelerator physics tools package for the OMC team at CERN."
 __url__ = "https://github.com/pylhc/omc3"
-__version__ = "0.6.6"
+__version__ = "0.7.0"
 __author__ = "pylhc"
 __author_email__ = "[email protected]"
 __license__ = "MIT"

omc3/amplitude_detuning_analysis.py

@@ -115,29 +115,27 @@
 
 
 """
-import os
-from collections import OrderedDict
 from datetime import timedelta
 from pathlib import Path
 from typing import List, Sequence, Tuple, Dict, Any, Union
 
 import numpy as np
 import tfs
 from generic_parser import DotDict
-from generic_parser.entrypoint_parser import EntryPointParameters, entrypoint, save_options_to_config
+from generic_parser.entrypoint_parser import EntryPointParameters, entrypoint
 from numpy.typing import ArrayLike
 from tfs.frame import TfsDataFrame
 
-from omc3.definitions import formats
 from omc3.definitions.constants import PLANES
 from omc3.tune_analysis import fitting_tools, kick_file_modifiers, timber_extract
-from omc3.tune_analysis.bbq_tools import OutlierFilterOpt, MinMaxFilterOpt
+from omc3.tune_analysis.bbq_tools import OutlierFilterOpt, MinMaxFilterOpt, FilterOpts
 from omc3.tune_analysis.constants import (
     get_bbq_col,
     get_bbq_out_name,
     get_kick_out_name,
-    get_mav_col,
-    get_timber_bbq_key, INPUT_KICK, INPUT_PREVIOUS, CORRECTED,
+    get_timber_bbq_key,
+    get_natq_err_col,
+    INPUT_KICK, INPUT_PREVIOUS, CORRECTED,
 )
 from omc3.tune_analysis.kick_file_modifiers import (
     read_timed_dataframe,
@@ -266,7 +264,8 @@ def analyse_with_bbq_corrections(opt: DotDict) -> Tuple[TfsDataFrame, TfsDataFra
 
     opt, filter_opt = _check_analyse_opt(opt)
     kick_df, bbq_df = get_kick_and_bbq_df(kick=opt.kick, bbq_in=opt.bbq_in,
-                                          beam=opt.beam, filter_opt=filter_opt, output=opt.output)
+                                          beam=opt.beam,
+                                          filter_opt=filter_opt)
 
     kick_plane = opt.plane
 
@@ -277,23 +276,34 @@ def analyse_with_bbq_corrections(opt: DotDict) -> Tuple[TfsDataFrame, TfsDataFra
             kick_df = double_action_analysis(kick_df, opt.detuning_order, corrected)
 
     if opt.output:
-        LOG.info(f"Writing kick data to file in directory '{opt.output.absolute()}'")
-        opt.output.mkdir(parents=True, exist_ok=True)
-        write_timed_dataframe(opt.output / get_kick_out_name(), kick_df)
-
+        _write_dataframes(opt.output, kick_df, bbq_df)
     return kick_df, bbq_df
 
 
 def get_kick_and_bbq_df(kick: Union[Path, str], bbq_in: Union[Path, str],
-                        beam: int = None, filter_opt = None,
-                        output: Path = None
+                        beam: int = None,
+                        filter_opt: FilterOpts = None,
                         ) -> Tuple[tfs.TfsDataFrame, tfs.TfsDataFrame]:
     """Load the input data."""
     bbq_df = None
     if bbq_in is not None and bbq_in == INPUT_PREVIOUS:
+        # NOTE: this is not the same as the "previous BBQ data" option in the GUI.
+        # That one just uses the previous bbq_ampdet.tfs file (loaded in the "else" below).
+        # The use-case for the INPUT_PREVIOUS option here is,
+        # that you can modify the kick_ampdet_xy file manually (e.g. removing kicks)
+        # and run the fitting on the new data again,
+        # without having to touch the whole BBQ stuff again (as the values are already in the file).
+        # Tips:
+        #  - Remove full columns to get rid of the whole kick
+        #  - Add NaNs into NATQ columns you want to ignore (in case you want to keep the other plane for this kick)
+        #  - Add NaNs to the ERRNATQ columns if you want to plot the point (w/o error bars) but not use it for fit
         LOG.debug("Getting data from previous ampdet kick file")
         kick_df = read_timed_dataframe(Path(kick) / get_kick_out_name())
         kick_df.headers = {k: v for k, v in kick_df.headers.items() if not k.startswith("ODR_")}
+
+        # redo the corrected columns, so you only need to add NaNs into the NATQ columns
+        LOG.debug("Adding corrected natural tunes and stdev to kick data")
+        kick_df = kick_file_modifiers.add_corrected_natural_tunes(kick_df)
     else:
         LOG.debug("Getting data from kick files")
         kick_df = read_two_kick_files_from_folder(kick)
@@ -307,15 +317,6 @@ def get_kick_and_bbq_df(kick: Union[Path, str], bbq_in: Union[Path, str],
             LOG.debug("Adding corrected natural tunes and stdev to kick data")
             kick_df = kick_file_modifiers.add_corrected_natural_tunes(kick_df)
 
-            if output:
-                LOG.info(f"Writing BBQ data to file in directory '{output.absolute()}'")
-                try:
-                    window = filter_opt.window
-                except AttributeError:
-                    window = filter_opt[0].window
-
-                x_interval = get_approx_bbq_interval(bbq_df, kick_df.index, window)
-                write_timed_dataframe(output / get_bbq_out_name(), bbq_df.loc[x_interval[0]: x_interval[1]])
     return kick_df, bbq_df
 
 
@@ -419,7 +420,7 @@ def get_approx_bbq_interval(
 # Private Functions ------------------------------------------------------------
 
 
-def _check_analyse_opt(opt: DotDict):
+def _check_analyse_opt(opt: DotDict) -> Tuple[DotDict, FilterOpts]:
     """Perform manual checks on opt-sturcture."""
     LOG.debug("Checking Options.")
 
@@ -542,6 +543,15 @@ def _get_ampdet_data_as_array(data: Dict[Any, AmpDetData], column: str) -> Array
     return np.vstack([getattr(d, column) for d in data.values()])
 
 
+def _write_dataframes(output: Path, kick_df: TfsDataFrame, bbq_df: TfsDataFrame):
+    LOG.info(f"Writing kick data to file in directory '{output.absolute()}'")
+    output.mkdir(parents=True, exist_ok=True)
+    write_timed_dataframe(output / get_kick_out_name(), kick_df)
+    if bbq_df is not None:
+        LOG.info(f"Writing BBQ data to file in directory '{output.absolute()}'")
+        write_timed_dataframe(output / get_bbq_out_name(), bbq_df)
+
+
 # Script Mode ##################################################################
 
 

omc3/harpy/handler.py

@@ -159,7 +159,7 @@ def _compute_headers(panda, date=None):
                 pass
             else:
                 headers[f"{prefix}Q{PLANE_TO_NUM[plane]}"] = np.mean(bpm_tunes)
-                headers[f"{prefix}Q{PLANE_TO_NUM[plane]}RMS"] = np.std(bpm_tunes)
+                headers[f"{prefix}Q{PLANE_TO_NUM[plane]}RMS"] = np.std(bpm_tunes)  # TODO: not really the RMS?
     if date:
         headers["TIME"] = date.strftime(formats.TIME)
     return headers

omc3/optics_measurements/constants.py

@@ -17,6 +17,7 @@
 ORBIT_NAME: str = "orbit_"
 KICK_NAME: str = "kick_"
 IP_NAME: str = "interaction_point_"
+CALIBRATION_FILE: str = "calibration_{plane}.out"
 
 # Column Names -----------------------------------------------------------------
 # Pre- and Suffixe
@@ -40,20 +41,28 @@
 DPP: str = "DPP"
 DPPAMP: str = "DPPAMP"
 AMPLITUDE: str = "AMP"
+NAT_AMPLITUDE: str = "NATAMP"
 PHASE: str = "PHASE"
 PHASE_ADV: str = "MU"
 F1001: str = "F1001"
 F1010: str = "F1010"
+NOISE: str = "NOISE"
+CLOSED_ORBIT: str = "CO"
 
 SECONDARY_AMPLITUDE_X: str = "AMP01_X"  # amplitude of secondary line in horizontal spectrum
 SECONDARY_AMPLITUDE_Y: str = "AMP10_Y"  # amplitude of secondary line in vertical spectrum
 SECONDARY_FREQUENCY_X: str = "PHASE01_X"  # frequency of secondary line in horizontal spectrum
 SECONDARY_FREQUENCY_Y: str = "PHASE10_Y"  # frequency of secondary line in vertical spectrum
 
+# Kick files
 TIME: str = "TIME"
 ACTION: str = "2J"
 SQRT_ACTION: str = "sqrt2J"
 
+# Calibration files
+CALIBRATION = "CALIBRATION"
+ERR_CALIBRATION = "ERROR_CALIBRATION"
 
 # Headers ----------------------------------------------------------------------
 RESCALE_FACTOR: str = "RescalingFactor"
+BPM_RESOLUTION: str = "BPMResolution"

omc3/optics_measurements/kick.py

@@ -20,7 +20,8 @@
                                                 NAT_TUNE, PEAK2PEAK,
                                                 RES,
                                                 RESCALE_FACTOR, RMS,
-                                                SQRT_ACTION, TIME, TUNE, S)
+                                                SQRT_ACTION, TIME, TUNE, S, NOISE, CLOSED_ORBIT)
+from omc3.utils.stats import weighted_mean, weighted_error
 
 
 def calculate(measure_input, input_files, scale, header_dict, plane):
@@ -37,7 +38,7 @@ def calculate(measure_input, input_files, scale, header_dict, plane):
         `TfsDataFrame` containing actions and their errors.
     """
     try:
-        kick_frame = _getkick(measure_input, input_files, plane)
+        kick_frame = _get_kick(measure_input, input_files, plane)
     except IndexError:  # occurs if either no x or no y files exist
         return pd.DataFrame
     kick_frame = _rescale_actions(kick_frame, scale, plane)
@@ -58,7 +59,7 @@ def _rescale_actions(df, scaling_factor, plane):
     return df
 
 
-def _getkick(measure_input, files, plane):
+def _get_kick(measure_input, files, plane):
     load_columns, calc_columns, column_types = _get_column_mapping(plane)
     kick_frame = pd.DataFrame(data=0.,
                               index=range(len(files[plane])),
@@ -71,22 +72,54 @@ def _getkick(measure_input, files, plane):
                 kick_frame.loc[i, col] = df[src]
 
         # calculate data from measurement
-        kick_frame.loc[i, calc_columns] = _gen_kick_calc(measure_input, df, plane)
+        kick_frame.loc[i, calc_columns] = _get_action(measure_input, df, plane)
     kick_frame = kick_frame.astype(column_types)
     return kick_frame
 
 
-def _gen_kick_calc(meas_input, lin, plane):
+def _get_action(meas_input, lin: pd.DataFrame, plane: str) -> np.ndarray:
     """
-    Takes either PK2PK/2 for kicker excitation or AMP for AC-dipole excitation
+    Calculates action (2J and sqrt(2J)) and its errors from BPM data in lin-df.
+    Takes either PK2PK/2 for kicker excitation or AMP for AC-dipole excitation,
+    as the amplitude of the oscillations for single kicks falls off over turns,
+    and hence the amplitude of the main line does not represent the initial kick,
+    whereas it is constant for the driven excitation.
+    Reminder: A = sqrt(2J \beta) .
+
+    TODO (jdilly 07.09.2022):
+          beta_phase instead of beta_model as stated below Eq. (11) in
+          PHYS. REV. ACCEL. BEAMS 23, 042801 (2020)
+
+    Returns:
+        sqrt(2J), error sqrt(2J), 2J, error 2J as  (1x4) array
     """
-    frame = pd.merge(_get_model_arc_betas(meas_input, plane), lin.loc[:, [f"{AMPLITUDE}{plane}", PEAK2PEAK]],
+    frame = pd.merge(_get_model_arc_betas(meas_input, plane), lin,
                      how='inner', left_index=True, right_index=True)
-    amps = (frame.loc[:, f"{AMPLITUDE}{plane}"].to_numpy() if meas_input.accelerator.excitation
-            else frame.loc[:, PEAK2PEAK].to_numpy() / 2.0)
-    meansqrt2j = amps / np.sqrt(frame.loc[:, f"{BETA}{plane}"].to_numpy())
-    mean2j = np.square(amps) / frame.loc[:, f"{BETA}{plane}"].to_numpy()
-    return np.array([np.mean(meansqrt2j), np.std(meansqrt2j), np.mean(mean2j), np.std(mean2j)])
+
+    if meas_input.accelerator.excitation:
+        amps = frame.loc[:, f"{AMPLITUDE}{plane}"].to_numpy()
+        err_amps = frame.loc[:, f"{ERR}{AMPLITUDE}{plane}"].to_numpy()
+    else:
+        amps = frame.loc[:, PEAK2PEAK].to_numpy() / 2.0
+        err_amps = frame.loc[:, f"{CLOSED_ORBIT}{RMS}"].to_numpy()
+
+    # sqrt(2J) ---
+    sqrt_beta = np.sqrt(frame.loc[:, f"{BETA}{plane}"].to_numpy())
+
+    actions_sqrt2j = amps / sqrt_beta
+    errors_sqrt2j = err_amps / sqrt_beta
+
+    mean_sqrt2j = weighted_mean(data=actions_sqrt2j, errors=errors_sqrt2j)
+    err_sqrt2j = weighted_error(data=actions_sqrt2j, errors=errors_sqrt2j)
+
+    # 2J ---
+    actions_2j = np.square(amps) / frame.loc[:, f"{BETA}{plane}"].to_numpy()
+    errors_2j = 2 * amps * err_amps / frame.loc[:, f"{BETA}{plane}"].to_numpy()
+
+    mean_2j = weighted_mean(data=actions_2j, errors=errors_2j)
+    err_2j = weighted_error(data=actions_2j, errors=errors_2j)
+
+    return np.array([mean_sqrt2j, err_sqrt2j, mean_2j, err_2j])
 
 
 def _get_model_arc_betas(measure_input, plane):

omc3/optics_measurements/measure_optics.py

@@ -10,6 +10,7 @@
 import sys
 from collections import OrderedDict
 from copy import deepcopy
+from typing import Dict, List
 
 import numpy as np
 import pandas as pd
@@ -21,7 +22,9 @@
                                       chromatic, dispersion, dpp, iforest,
                                       interaction_point, kick, phase, rdt,
                                       tune, crdt, coupling)
-from omc3.optics_measurements.constants import (CHROM_BETA_NAME, ERR, EXT)
+from omc3.optics_measurements.constants import (
+    CHROM_BETA_NAME, ERR, EXT, AMPLITUDE, ERR_CALIBRATION, CALIBRATION, CALIBRATION_FILE, NAME, BPM_RESOLUTION
+)
 from omc3.utils import iotools, logging_tools
 
 LOGGER = logging_tools.get_logger(__name__, level_console=logging_tools.INFO)
@@ -47,7 +50,7 @@ def measure_optics(input_files, measure_input):
     invariants = {}
     phase_dict = {}
     for plane in PLANES:
-        phase_dict[plane], out_dfs  = phase.calculate(measure_input, input_files, tune_dict, plane)
+        phase_dict[plane], out_dfs = phase.calculate(measure_input, input_files, tune_dict, plane)
         phase.write(out_dfs, [common_header]*4, measure_input.outputdir, plane)
         phase.write_special(measure_input, phase_dict[plane]['free'], tune_dict[plane]["QF"], plane)
         if measure_input.only_coupling:
@@ -150,7 +153,8 @@ def __init__(self, files_to_analyse, optics_opt):
     def _repair_backwards_compatible_frame(df, plane):
         """
         Multiplies unscaled amplitudes by 2 to get from complex amplitudes to the real ones.
-        This is for backwards compatibility with Drive.
+        This is for backwards compatibility with Drive,
+        i.e. harpy has this
         """
         df[f"AMP{plane}"] = df.loc[:, f"AMP{plane}"].to_numpy() * 2
         if f"NATAMP{plane}" in df.columns:
@@ -220,16 +224,34 @@ def bpms(self, plane=None, dpp_value=None):
             indices[0] = indices[0].intersection(ind)
         return indices[0]
 
-    def calibrate(self, calibs):
+    def calibrate(self, calibs: Dict[str, pd.DataFrame]):
+        """
+        Use calibration data to rescale amplitude and amplitude error.
+
+        Args:
+            calibs (Dict): Plane-Dictionary with DataFrames of calibration data.
+
+        """
         if calibs is None:
             return
+
         for plane in PLANES:
+            bpm_resolution = calibs[plane].headers.get(BPM_RESOLUTION, 1e-4)  # TODO: Default of 0.1 mm is LHC specific
             for i in range(len(self[plane])):
-                data = pd.merge(self[plane][i].loc[:, ["AMP" + plane]], calibs[plane], how='left',
-                                left_index=True, right_index=True).fillna(
-                    value={"CALIBRATION": 1., "ERROR_CALIBRATION": 0.})
-                self[plane][i][f"AMP{plane}"] = self[plane][i].loc[:, f"AMP{plane}"] * data.loc[:, "CALIBRATION"]
-                self[plane][i][f"{ERR}AMP{plane}"] = data.loc[:, "ERROR_CALIBRATION"]  # TODO
+                # Merge all measurement BPMs into calibration data (only few BPMs),
+                # fill missing values with a scaling of 1 and estimated error of 0.1mm (emaclean estimate)
+                data = pd.merge(self[plane][i].loc[:, [f"{AMPLITUDE}{plane}"]], calibs[plane],
+                                how='left', left_index=True, right_index=True).fillna(
+                    value={CALIBRATION: 1.})  # ERR_CALIBRATION is relative, NaN filled with absolute value below
+
+                # Scale amplitude with the calibration
+                self[plane][i][f"AMP{plane}"] = self[plane][i].loc[:, f"AMP{plane}"] * data.loc[:, CALIBRATION]
+
+                # Sum Amplitude Error (absolute) and Calibration Error (relative)
+                self[plane][i][f"{ERR}{AMPLITUDE}{plane}"] = np.sqrt(
+                    self[plane][i][f"{ERR}{AMPLITUDE}{plane}"]**2 +
+                    ((self[plane][i][f"{AMPLITUDE}{plane}"] * data.loc[:, ERR_CALIBRATION]).fillna(bpm_resolution))**2
+                )
 
     @ staticmethod
     def get_columns(frame, column):
@@ -268,7 +290,7 @@ def copy_calibration_files(outputdir, calibrationdir):
         return None
     calibs = {}
     for plane in PLANES:
-        cal_file = f"calibration_{plane.lower()}.out"
+        cal_file = CALIBRATION_FILE.format(plane=plane.lower())
         iotools.copy_item(os.path.join(calibrationdir, cal_file), os.path.join(outputdir, cal_file))
-        calibs[plane] = tfs.read(os.path.join(outputdir, cal_file)).set_index("NAME")
+        calibs[plane] = tfs.read(os.path.join(outputdir, cal_file)).set_index(NAME)
     return calibs

omc3/scripts/betabeatsrc_output_converter.py

@@ -1,6 +1,6 @@
 """
 BetaBeat.src Output Converter
--------------
+-----------------------------
 
 Script to convert most important output files produced by ``BetaBeat.src`` / ``GetLLM`` into the standard
 format used in ``omc3`` to allow straight forward comparison of the two.