diff --git a/common/src/aaredaqlib/find_xtal.py b/common/src/aaredaqlib/find_xtal.py new file mode 100644 index 00000000..1eb493d9 --- /dev/null +++ b/common/src/aaredaqlib/find_xtal.py @@ -0,0 +1,284 @@ +from typing import List, Optional, Callable + +import numpy as np +from scipy import ndimage + +from aaredaqlib.models import CrystalSize +from aaredaqlib.raster_grid import RasterGridRequest + + +def identify_crystal_raster(result, r: RasterGridRequest): + images = result.images + if images and any(img.spots > 0 for img in images): + + indexed_images = [img for img in images if img.index and img.spots_low_res > 4 and img.bkg > 4.5] + + if indexed_images: + # indexed_images = [img for img in indexed_images if img.spots_indexed > 10] + images = indexed_images + + filtered_images = [img for img in images + if img.spots_ice is not None and img.spots_low_res > 4 and ( + img.spots_ice / img.spots_low_res) < 5.0 + and (img.spots_ice / img.spots_low_res) != 1] + + if indexed_images: + print(f"Find image by maximum number of spots indexed") + max_image = max(images, key=lambda img: img.spots_indexed) + max_spots = max_image.spots_indexed + max_images = [img for img in images if img.spots_indexed == max_spots] + max_image = max_images[len(max_images) // 2] + + print(f"Image with maximum spots_low_res: {max_image}") + print(f"Maximum spots_indexed value: {max_image.spots_indexed}") + print(f"Maximum spots_low_res value: {max_image.spots_low_res}") + else: + print(f"Find image by maximum number of low resolution spots") + max_image = max(images, key=lambda img: img.spots_low_res) + print(f"Image with maximum spots_low_res: {max_image}") + print(f"Maximum spots_low_res value: {max_image.spots_low_res}") + + grid_mm_x = max_image.nx * r.grid_size_mm.x + grid_mm_y = max_image.ny * r.grid_size_mm.y + + print(f"Grid coordinates in mm: ({grid_mm_x}, {grid_mm_y})") + return grid_mm_x, grid_mm_y + + else: + return None, None + +def rebuild_array_from_scan_results(scan_results: List, + value_field: str, + array_shape: Optional[tuple] = None, + nx_field: str = 'nx', + ny_field: str = 'ny', + default_value: float = 0.0, + threshold: Optional[float] = None, + condition_func: Optional[Callable] = None, + apply_filter_before: bool = True + ) -> np.ndarray: + # Extract coordinates and values + positions = [] + values = [] + + for result in scan_results: + nx = getattr(result, nx_field) + ny = getattr(result, ny_field) + value = getattr(result, value_field) + + # Skip if coordinates are None + if nx is None or ny is None: + continue + + positions.append((int(nx), int(ny))) # Note: (row, col) = (ny, nx) + if not value: + value = 0.0 + values.append(float(value)) + + if not positions: + raise ValueError("No valid positions found in scan results") + + # Determine array shape + if array_shape is None: + max_row = max(pos[0] for pos in positions) + max_col = max(pos[1] for pos in positions) + array_shape = (max_row + 1, max_col + 1) + + # Initialize array with default values + result_array = np.full(array_shape, default_value, dtype=float) + + # Apply pre-filtering if requested + if apply_filter_before: + filtered_data = [] + for pos, val in zip(positions, values): + keep_value = True + + # Apply threshold filter + if threshold is not None and val < threshold: + keep_value = False + + # Apply custom condition + if condition_func is not None and not condition_func(val): + keep_value = False + + if keep_value: + filtered_data.append((pos, val)) + else: + filtered_data.append((pos, 0.0)) + + # Fill array with filtered values + for pos, val in filtered_data: + if 0 <= pos[0] < array_shape[0] and 0 <= pos[1] < array_shape[1]: + result_array[pos[0], pos[1]] = val + else: + # Fill array first, then apply filters + for pos, val in zip(positions, values): + if 0 <= pos[0] < array_shape[0] and 0 <= pos[1] < array_shape[1]: + result_array[pos[0], pos[1]] = val + + # Apply post-filtering + if threshold is not None: + result_array[result_array < threshold] = 0.0 + + if condition_func is not None: + mask = np.vectorize(condition_func)(result_array) + result_array[~mask] = 0.0 + + return result_array + +def create_quality_filtered_array(scan_results: List, + value_field: str, + min_spots: Optional[int] = None, + min_efficiency: Optional[float] = 1.0, + min_background: Optional[float] = None, + exclude_ice: Optional[bool] = True, + min_low_res_spots: Optional[float] = 10.0, + **kwargs + ) -> np.ndarray: + """ + Create array with comprehensive quality filtering + """ + + def quality_condition(result, min_bkg, min_spots, min_efficiency, min_res_spots = 10.0): + if exclude_ice and (result.spots_ice / max(result.spots_low_res, 1.0)) == 1.0: + # print(f"all ice for {result.number}") + return False + if result.spots_low_res < min_res_spots: + return False + if exclude_ice and result.spots_ice > result.spots * 0.8: # More than 50% ice + # print(f"more than 80% ice for {result.number}") + return False + if result.index: + # print(f"index is True for {result.number}") + return True + if result.spots < min_spots: + # print(f"{result.spots} is less than {min_spots} for {result.number}") + return False + if result.spots_low_res < min_background: + return False + + if result.efficiency < min_efficiency: + return False + + return True + + # Filter results first + filtered_results = [] + + if min_spots is None: + min_spots = min((result.spots for result in scan_results if result.spots is not None), default=1) + if min_background is None: + min_background = min((result.bkg for result in scan_results if result.bkg is not None), default=1) + if min_efficiency is None: + min_efficiency = 1.0 + + for result in scan_results: + if result.nx is not None and result.ny is not None: + if quality_condition(result, min_bkg=min_background, min_spots=min_spots, + min_efficiency=min_efficiency, min_res_spots=min_low_res_spots): + filtered_results.append(result) + else: + # Create a copy with zero value for filtered positions + import copy + zero_result = copy.copy(result) + setattr(zero_result, value_field, 0) + filtered_results.append(zero_result) + + return rebuild_array_from_scan_results(filtered_results, value_field, **kwargs) + + +def get_xtal_size(crystal_size, result_array, r:RasterGridRequest): + # Optional: get bounding box of the largest object + try: + labeled_array, num_objects = ndimage.label(result_array) + areas = ndimage.sum(np.ones_like(result_array, dtype=np.int32), labeled_array, + index=range(1, num_objects + 1)) + largest_idx = int(np.argmax(areas)) + 1 # +1 because labels start at 1 + largest_area = int(areas[largest_idx - 1]) + print(f"Largest object label: {largest_idx}, area (px): {largest_area}") + + object_mask = labeled_array == largest_idx + rows = np.any(object_mask, axis=1) + cols = np.any(object_mask, axis=0) + row_min, row_max = np.where(rows)[0][[0, -1]] + col_min, col_max = np.where(cols)[0][[0, -1]] + print(f"Largest bbox: width={col_max - col_min}, height={row_max - row_min}") + print( + f"Largest bbox: width={(col_max - col_min) * r.grid_size_mm.x}, y={(row_max - row_min) * r.grid_size_mm.y}") + if r.n_x == 1: + crystal_size = CrystalSize(x=crystal_size.x, y=crystal_size.y, + z=(col_max - col_min) * r.grid_size_mm.y * 1000) + + else: + crystal_size = CrystalSize(x=(row_max - row_min) * r.grid_size_mm.x * 1000, + y=(col_max - col_min) * r.grid_size_mm.y * 1000, + z=crystal_size.z) + except ValueError as e: + print(f"error calculating xtal size: {e}") + crystal_size = CrystalSize(x=0,y=0,z=0) + + return crystal_size + + +def get_best_b_factor(result_list: List): + if not result_list: + return None + best_b_factor = min((img for img in result_list if img.b is not None), + key=lambda img: img.b, + default=None) + if best_b_factor is None: + return None + print(f"Best b: {best_b_factor.b}") + return best_b_factor.b + +def get_best_res(result_list: List): + if not result_list: + return None + best_res = min((img for img in result_list if img.res is not None), + key=lambda img: img.res, + default=None) + print(f"Best res: {best_res.res}") + return best_res.res + +def com_nan_check(com): + if np.isnan(com[0]) or np.isnan(com[1]): + return False + else: + return True + +def get_result_list_from_com(images, com): + if not com_nan_check(com): + return None + cy, cx = com[::-1] + start_x, end_x = round(cx - 1), round(cx + 1) + start_y, end_y = round(cy - 1), round(cy + 1) + print(f"range x {start_x} {end_x}, y {start_y} {end_y}") + result_list = [img for img in images + if start_x <= img.nx <= end_x and start_y <= img.ny <= end_y] + return result_list + +def get_com_image_number(com, images): + for image in images: + if image.nx == round(com[0]) and image.ny == round(com[1]): + print(f"com found for image: {image.number}") + return + +def get_grid_mm_from_com(com, r:RasterGridRequest): + if not com_nan_check(com): + grid_mm_x = None + grid_mm_y = None + else: + grid_mm_x = com[0] * r.grid_size_mm.x + grid_mm_y = com[1] * r.grid_size_mm.y + if r.n_x == 1: + grid_mm_x += (0.5 * r.grid_size_mm.x) + return grid_mm_x, grid_mm_y + + +def raster_centre_of_mass(result_array, r:RasterGridRequest): + print('horizontal scan') + com = ndimage.center_of_mass(result_array) + print(f"Center of mass: {com}") + grid_mm_x, grid_mm_y = get_grid_mm_from_com(com, r) + return grid_mm_x, grid_mm_y, com + diff --git a/common/src/aaredaqlib/models.py b/common/src/aaredaqlib/models.py index 3beb5964..bf9b55fe 100644 --- a/common/src/aaredaqlib/models.py +++ b/common/src/aaredaqlib/models.py @@ -403,6 +403,13 @@ class SampleShortInfo(BaseModel): class SampleShortInfoList(BaseModel): s: List[SampleShortInfo] +class SimpleScanParameters(BaseModel): + dtz: float = 120 + exp_time_s: float = 0.02 + start_omega_deg:float = 0 + incr_omega_deg: float = 0.2 + steps: int = 900 + transmission: float = 1.0 #class ReferencePuckInfo(BaseModel): # test_sample = SampleShortInfo(db_id=-1, puck_name="test_puck", dewar_name="test_dewar", sample_name="test_sample", # pin=11, location=DewarAddress(segment="X", pos=1)) diff --git a/daq/src/aaredaq/best.pt b/daq/src/aaredaq/best.pt index 8ad56d3d..ed71546a 100644 Binary files a/daq/src/aaredaq/best.pt and b/daq/src/aaredaq/best.pt differ diff --git a/daq/src/aaredaq/daq.py b/daq/src/aaredaq/daq.py index 8bb68359..3270e9e0 100644 --- a/daq/src/aaredaq/daq.py +++ b/daq/src/aaredaq/daq.py @@ -6,6 +6,7 @@ from datetime import datetime from math import ceil from typing import List, Tuple, Optional, Callable, Dict import secrets +import os import cv2 import numpy as np @@ -23,12 +24,14 @@ from aaredaqlib.beamline import MXBeamline from aaredaqlib.coordinate import Coordinate, SmargonCoordinate from aaredaqlib.diffraction_geometry import DiffractionGeometry from aaredaqlib.logger_config import setup_logger +from aaredaqlib.find_xtal import raster_centre_of_mass, create_quality_filtered_array, identify_crystal_raster, \ + get_result_list_from_com, get_best_b_factor, get_best_res, get_xtal_size from aaredaqlib.models import ( SampleShortInfo, PuckLoadedInfo, SampleShortInfoList, DAQStatusModel, BeamlineStatus, SessionStatus, SampleCameraSettings, AutofocusSettings, BoundingBoxModel, - ZoomModeEnum, CrystalSize) + ZoomModeEnum, CrystalSize, SimpleScanParameters) from aaredaqlib.raster_grid import RasterGridRequest, CompletedRasterGrid, CompletedRasterGridElem from aaredaqlib.rotation_scan import RotationScanRequest, CompletedRotationScan from aaredaqlib.sample_geometry import SampleGeometryModel @@ -64,6 +67,7 @@ class AareDAQ: self.crystal_size = CrystalSize(x=0,y=0,z=0) self.last_best_b_factor = None self.last_best_res = None + self.auto_params = SimpleScanParameters() @property def state(self) -> BeamlineStateEnum: @@ -88,6 +92,14 @@ class AareDAQ: self.last_time = end - start + @property + def smart_params(self) -> SimpleScanParameters: + return self.auto_params + + @smart_params.setter + def smart_params(self, params: SimpleScanParameters): + self.auto_params = params + @property def omega(self) -> float: return self.__devs.aerotech.omega @@ -369,202 +381,22 @@ class AareDAQ: def list_loaded_pucks(self) -> List[PuckLoadedInfo]: return self.__devs.tell.get_detected_pucks() - def identify_crystal_raster(self, result, r: RasterGridRequest): - images = result.images - if images and any(img.spots > 0 for img in images): - - indexed_images = [img for img in images if img.index and img.spots_low_res > 4 and img.bkg > 4.5] - - if indexed_images: - #indexed_images = [img for img in indexed_images if img.spots_indexed > 10] - images = indexed_images - - filtered_images = [img for img in images - if img.spots_ice is not None and img.spots_low_res > 4 and ( - img.spots_ice / img.spots_low_res) < 5.0 - and (img.spots_ice / img.spots_low_res) != 1] - - if indexed_images: - print(f"Find image by maximum number of spots indexed") - max_image = max(images, key=lambda img: img.spots_indexed) - max_spots = max_image.spots_indexed - max_images = [img for img in images if img.spots_indexed == max_spots] - max_image = max_images[len(max_images) // 2] - - print(f"Image with maximum spots_low_res: {max_image}") - print(f"Maximum spots_indexed value: {max_image.spots_indexed}") - print(f"Maximum spots_low_res value: {max_image.spots_low_res}") - else: - print(f"Find image by maximum number of low resolution spots") - max_image = max(images, key=lambda img: img.spots_low_res) - print(f"Image with maximum spots_low_res: {max_image}") - print(f"Maximum spots_low_res value: {max_image.spots_low_res}") - - grid_mm_x = max_image.nx * r.grid_size_mm.x - grid_mm_y = max_image.ny * r.grid_size_mm.y - - print(f"Grid coordinates in mm: ({grid_mm_x}, {grid_mm_y})") - delta_mm = self.sample_geometry.smargon_nudge(Coordinate(x=grid_mm_x, y=grid_mm_y)) - return delta_mm - - else: - return None - - def rebuild_array_from_scan_results(self, - scan_results: List, - value_field: str, - array_shape: Optional[tuple] = None, - nx_field: str = 'nx', - ny_field: str = 'ny', - default_value: float = 0.0, - threshold: Optional[float] = None, - condition_func: Optional[Callable] = None, - apply_filter_before: bool = True - ) -> np.ndarray: - - # Extract coordinates and values - positions = [] - values = [] - - for result in scan_results: - nx = getattr(result, nx_field) - ny = getattr(result, ny_field) - value = getattr(result, value_field) - - # Skip if coordinates are None - if nx is None or ny is None: - continue - - positions.append((int(nx), int(ny))) # Note: (row, col) = (ny, nx) - if not value: - value = 0.0 - values.append(float(value)) - - if not positions: - raise ValueError("No valid positions found in scan results") - - # Determine array shape - if array_shape is None: - max_row = max(pos[0] for pos in positions) - max_col = max(pos[1] for pos in positions) - array_shape = (max_row + 1, max_col + 1) - - # Initialize array with default values - result_array = np.full(array_shape, default_value, dtype=float) - - # Apply pre-filtering if requested - if apply_filter_before: - filtered_data = [] - for pos, val in zip(positions, values): - keep_value = True - - # Apply threshold filter - if threshold is not None and val < threshold: - keep_value = False - - # Apply custom condition - if condition_func is not None and not condition_func(val): - keep_value = False - - if keep_value: - filtered_data.append((pos, val)) - else: - filtered_data.append((pos, 0.0)) - - # Fill array with filtered values - for pos, val in filtered_data: - if 0 <= pos[0] < array_shape[0] and 0 <= pos[1] < array_shape[1]: - result_array[pos[0], pos[1]] = val - else: - # Fill array first, then apply filters - for pos, val in zip(positions, values): - if 0 <= pos[0] < array_shape[0] and 0 <= pos[1] < array_shape[1]: - result_array[pos[0], pos[1]] = val - - # Apply post-filtering - if threshold is not None: - result_array[result_array < threshold] = 0.0 - - if condition_func is not None: - mask = np.vectorize(condition_func)(result_array) - result_array[~mask] = 0.0 - - return result_array - - def create_quality_filtered_array(self, - scan_results: List, - value_field: str, - min_spots: Optional[int] = None, - min_efficiency: Optional[float] = 1.0, - min_background: Optional[float] = None, - exclude_ice: Optional[bool] = True, - **kwargs - ) -> np.ndarray: - - """ - Create array with comprehensive quality filtering - """ - - def quality_condition(result, min_bkg, min_spots, min_efficiency): - if exclude_ice and (result.spots_ice / max(result.spots_low_res, 1.0)) == 1.0: - #print(f"all ice for {result.number}") - return False - # if (result.spots_ice / result.spots_low_res) > 5.0: - # return False - if exclude_ice and result.spots_ice > result.spots * 0.8: # More than 50% ice - #print(f"more than 80% ice for {result.number}") - return False - if result.index: - #print(f"index is True for {result.number}") - return True - if result.spots < min_spots: - # print(f"{result.spots} is less than {min_spots} for {result.number}") - return False - if result.spots_low_res < min_background: - return False - - if result.efficiency < min_efficiency: - return False - - return True - - # Filter results first - filtered_results = [] - - if min_spots is None: - min_spots = min((result.spots for result in scan_results if result.spots is not None), default=1) - if min_background is None: - min_background = min((result.bkg for result in scan_results if result.bkg is not None), default=1) - if min_efficiency is None: - min_efficiency = 1.0 - - for result in scan_results: - if result.nx is not None and result.ny is not None: - if quality_condition(result, min_bkg=min_background, min_spots=min_spots, - min_efficiency=min_efficiency): - filtered_results.append(result) - else: - # Create a copy with zero value for filtered positions - import copy - zero_result = copy.copy(result) - setattr(zero_result, value_field, 0) - filtered_results.append(zero_result) - - return self.rebuild_array_from_scan_results(filtered_results, value_field, **kwargs) - def __auto_center(self, grid: RasterGridRequest) -> CompletedRasterGrid | None: sample = self.sample + if sample is None: raise Exception("Sample must be mounted to auto center") old_prefix = grid.file_prefix geom = self.sample_geometry r = self.__ml_bounding_box(sample.db_id, f"ml_{geom.omega_deg:.2f}deg") + if r is None: self.__devs.aerotech.move(geom.omega_deg + 90.0, wait=True) time.sleep(0.2) r = self.__ml_bounding_box(sample.db_id, f"ml_{geom.omega_deg + 90.0:.2f}deg") + if r is not None: geom = self.sample_geometry grid.smargon = r.smargon @@ -590,132 +422,6 @@ class AareDAQ: else: return None - def raster_centre_of_mass(self, images, r:RasterGridRequest, result): - # if any(img.index for img in images): - # print("COM by indexed spots") - # result_array= self.create_quality_filtered_array(images, 'spots_indexed', min_spots=None, - # min_efficiency=1.0, min_background=None) - # else: - print("COM by low res spots") - result_array = self.create_quality_filtered_array(images, 'spots_low_res', min_spots=None, - min_efficiency=1.0, min_background=None) - - print('horizontal scan') - com = ndimage.center_of_mass(result_array) - - print(f"Center of mass: {com}") - - try: - labeled_array, num_objects = ndimage.label(result_array) - areas = ndimage.sum(np.ones_like(result_array, dtype=np.int32), labeled_array, - index=range(1, num_objects + 1)) - largest_idx = int(np.argmax(areas)) + 1 # +1 because labels start at 1 - largest_area = int(areas[largest_idx - 1]) - print(f"Largest object label: {largest_idx}, area (px): {largest_area}") - - # Optional: get bounding box of largest object - object_mask = labeled_array == largest_idx - rows = np.any(object_mask, axis=1) - cols = np.any(object_mask, axis=0) - row_min, row_max = np.where(rows)[0][[0, -1]] - col_min, col_max = np.where(cols)[0][[0, -1]] - print(f"Largest bbox: width={col_max - col_min}, height={row_max - row_min}") - print( - f"Largest bbox: width={(col_max - col_min) * r.grid_size_mm.x}, y={(row_max - row_min) * r.grid_size_mm.y}") - - if r.n_x == 1: - - if self.crystal_size is None: - self.crystal_size = CrystalSize(x=0,y=0,z=0) - - crystal_size = self.crystal_size - - crystal_size = CrystalSize(x=crystal_size.x, y=crystal_size.y, - z=(col_max - col_min) * r.grid_size_mm.y * 1000) - - else: - crystal_size = CrystalSize(x=(row_max - row_min) * r.grid_size_mm.x * 1000, - y=(col_max - col_min) * r.grid_size_mm.y * 1000, - z=0) - except Exception as e: - crystal_size = CrystalSize(x=0, y=0, z=0) - - - self.crystal_size = crystal_size - - if r.n_x == 1 and (np.isnan(com[1]) or np.isnan(com[0])): - print('vertical scan') - try: - max_image = max(images, key=lambda img: img.spots_low_res) - com = (0, max_image.ny) - except: - print("no spots") - - if np.isnan(com[1]) or np.isnan(com[0]): - print("Center of mass is nan") - com = None - grid_mm_x = None - grid_mm_y = None - best_res = None - best_b_factor = None - - else: - grid_mm_x = com[0] * r.grid_size_mm.x - grid_mm_y = com[1] * r.grid_size_mm.y - if r.n_x == 1: - grid_mm_x += (0.5 * r.grid_size_mm.x) - - cx, cy = com[::-1] # com=(y, x) -> (x, y), rounded once - start_x, end_x = round(cx - 1), round(cx + 1) - start_y, end_y = round(cy - 1), round(cy + 1) - print(f"range x {start_x} {end_x}, y {start_y} {end_y}") - # Collect images in the 3x3 neighborhood around the center - res_list = [img for img in images - if start_x <= img.nx <= end_x and start_y <= img.ny <= end_y] - print(res_list) - # Best by res, skipping None - best_res = min((img for img in res_list if img.res is not None), - key=lambda img: img.res, - default=None) - print(f"Best res: {best_res}") - best_b_factor = min((img for img in res_list if img.b is not None), - key=lambda img: img.b, - default=None) - print(f"Best b: {best_b_factor}") - - for image in images: - if image.nx == round(com[0]) and image.ny == round(com[1]): - print(f"com found for image: {image.number}") - try: - if best_res is not None and best_res.res is not None: - # store the numeric resolution on the config/session so it appears in status - print(f"best res: {best_res.res}") - self.last_best_res = float(best_res.res) - else: - print(f'res is None') - self.last_best_res = None - except Exception as e: - print(f'error with last_best_res {e}') - self.last_best_res = None - try: - if best_b_factor is not None and best_b_factor.b is not None: - print(f"best res: {best_res.b}") - self.last_best_b_factor = float(best_b_factor.b) - else: - print(f'last_best_b_factor is None') - self.last_best_b_factor = None - except Exception as e: - print(f'error with last_best_b_factor {e}') - self.last_best_b_factor = None - - if com is not None and grid_mm_x is not None and grid_mm_y is not None: - new_delta_mm = self.sample_geometry.smargon_nudge(Coordinate(x=grid_mm_x, y=grid_mm_y)) - print(f"new delta mm: {new_delta_mm}, new grid x: {grid_mm_x}, new grid y: {grid_mm_y}") - else: - print(f"using old method as COM is none or nan") - new_delta_mm = self.identify_crystal_raster(result, r) - - return new_delta_mm, best_res def __raster(self, r: RasterGridRequest) -> CompletedRasterGridElem: max_time = r.exp_time_s * r.n_y * r.n_x + 60 @@ -774,7 +480,28 @@ class AareDAQ: with open(filename, 'w') as f: json.dump(output_data, f, indent=2) print('before centre_of_mass') - new_delta_mm, best_res = self.raster_centre_of_mass(images, r, result) + print("COM by low res spots") + + + result_array = create_quality_filtered_array(images, 'spots_low_res', min_spots=None, + min_efficiency=1.0, min_background=None, min_low_res_spots=10.0) + self.crystal_size = get_xtal_size(self.crystal_size, result_array, r) + grid_mm_x, grid_mm_y, com = raster_centre_of_mass(result_array, r) + + if grid_mm_x is None or grid_mm_y is None: + print(f"using old method as COM is none or nan") + grid_mm_x, grid_mm_y = identify_crystal_raster(result, r) + + if r.n_x == 1: + new_delta_mm = self.sample_geometry.smargon_nudge(Coordinate(x=0, y=grid_mm_y)) + else: + result_list = get_result_list_from_com(images, com) + print(result_list) + self.last_best_b_factor = get_best_b_factor(result_list) + self.last_best_res = get_best_res(result_list) + print("b_factor: ", self.last_best_b_factor, " best_res: ", self.last_best_res) + new_delta_mm = self.sample_geometry.smargon_nudge(Coordinate(x=grid_mm_x, y=grid_mm_y)) + print('after centre_of_mass') if new_delta_mm is not None: print(f'{time.ctime()}, moving SMARGON to target new delta mm {r.smargon.sh_mm + new_delta_mm} mm') @@ -1490,7 +1217,6 @@ class AareDAQ: def measure(self, sample: SampleShortInfo) -> float: start = time.perf_counter() - formatted_date = datetime.now().strftime('%Y%m%d') sample_prefix = "{}/{}/{:02d}/{}".format( formatted_date, @@ -1533,6 +1259,8 @@ class AareDAQ: grid_size_mm=Coordinate(x=geom.beam_size_mm.x * 0.5, y=geom.beam_size_mm.y * 0.5), omega_deg=0 )): + params = self.smart_params + print(params) self.__cfg.zoom_mode = ZoomModeEnum.User self.__devs.samcam_settings = self.__cfg.zoom_settings.get_camera_settings(self.zoom) self.__rotation( RotationScanRequest(start_omega_deg=0, @@ -1542,13 +1270,14 @@ class AareDAQ: incr_omega_deg=0.2, steps=900,)) else: + print("auto center failed") self.__aare.axc_failed(sample) self.zoom = 1 self.__cfg.zoom_mode = ZoomModeEnum.User self.__devs.samcam_settings = self.__cfg.zoom_settings.get_camera_settings(self.zoom) self.__cfg.state_busy = False - except Exception: - logger.exception("Exception in measure") + except Exception as e: + print(f"Error in measure: {e}") self.__aare.sample_failed(sample) self.__cfg.state_busy = False end = time.perf_counter() diff --git a/daq/src/aaredaq/server.py b/daq/src/aaredaq/server.py index 559e3667..2461d926 100644 --- a/daq/src/aaredaq/server.py +++ b/daq/src/aaredaq/server.py @@ -8,7 +8,7 @@ import uvicorn from aaredaqlib.coordinate import SmargonCoordinate, Coordinate from aaredaqlib.models import SampleShortInfo, DAQStatusModel, BeamlineStateEnum, BeamlineSettingsModel, \ SampleShortInfoList, SessionStatus, SampleCameraSettings, AutofocusSettings, TokenData, \ - CryojetSettingsModel + CryojetSettingsModel, SimpleScanParameters from aaredaqlib.raster_grid import RasterGridRequest, CompletedRasterGrid from aaredaqlib.rotation_scan import RotationScanRequest, CompletedRotationScan from aaredaqlib.sample_geometry import SampleGeometryModel @@ -345,6 +345,12 @@ async def auto(s: SampleShortInfo, token: str = Depends(oauth2_scheme)): runtime = daq.measure(s) return f"{runtime:0.3f}" +@app.post("/scan/smart_params") +async def set_smart_params(p: SimpleScanParameters, token: str = Depends(oauth2_scheme)) -> str: + auth.check_jwt_rw(cfg, auth.parse_token(token)) + daq.smart_params = p + return "OK" + @app.post("/scan/cancel") async def cancel(token: str = Depends(oauth2_scheme)): auth.check_jwt_rw(cfg, auth.parse_token(token)) diff --git a/daq/src/aaredaq/workflows.py b/daq/src/aaredaq/workflows.py index 5c49b82d..fe075ded 100644 --- a/daq/src/aaredaq/workflows.py +++ b/daq/src/aaredaq/workflows.py @@ -143,6 +143,7 @@ def dc2rse(devs: BeamlineDevices, cfg: BeamlineConfig): def se2sa(devs: BeamlineDevices, cfg: BeamlineConfig): hub = cfg.settings cfg.zoom_mode = ZoomModeEnum.User + devs.detector_cover.put(1) devs.samcam_settings = cfg.zoom_settings.get_camera_settings(devs.zoom) hub_cryo = cfg.cryojet_settings _cryo = hub_cryo.cryojet_in_use @@ -152,7 +153,7 @@ def se2sa(devs: BeamlineDevices, cfg: BeamlineConfig): devs.aerotech.unlock() devs.aerotech.set_direct_mode() print(time.ctime(), " moving reflector to up position") - devs.reflector_up = True + #devs.reflector_up = True print(time.ctime(), " moving beamstop to up position") devs.beamstop_stage_up = True print(time.ctime(), " setting lamp to 2.5") @@ -167,8 +168,6 @@ def se2sa(devs: BeamlineDevices, cfg: BeamlineConfig): devs.dtz.move(cfg.dtz, wait=False) def rse2sa(devs: BeamlineDevices, cfg: BeamlineConfig): - devs.detector_cover.put(2) # Open detector cover - print(time.ctime(), " moving DETECTOR COVER to OPEN position") if devs.tell.is_in_park(): print(time.ctime(), " moving TELL to COLD") devs.tell.move_cold(wait=True) diff --git a/daq/src/mxlibs3/jfjoch.py b/daq/src/mxlibs3/jfjoch.py index 68a82e9c..68c53b6e 100644 --- a/daq/src/mxlibs3/jfjoch.py +++ b/daq/src/mxlibs3/jfjoch.py @@ -126,8 +126,11 @@ class JFJochWrapper: ring_current_mA=s.bl.ring_current_mA, sample_temperature_K=s.bl.cryojet_K, total_flux=s.bl.flux_ph_s, - space_group_number=1 - #TODO space_group_number and unit_cell??? + space_group_number=1,#TODO space_group_number and unit_cell??? + poni_rot1_rad = s.diffraction.poni_rot1_rad, + poni_rot2_rad = s.diffraction.poni_rot2_rad, + max_spot_count = 1000, + detect_ice_rings = True ) self.__api.start_post(dataset_settings=dataset_settings) diff --git a/gui/src/aaregui/main_window.py b/gui/src/aaregui/main_window.py index 41cfdee5..96b501fc 100644 --- a/gui/src/aaregui/main_window.py +++ b/gui/src/aaregui/main_window.py @@ -68,8 +68,8 @@ class MainWindow(QMainWindow): beam_center_pxl=(750, 750), detector_description="PILATUS 4", detector_serial_number="1", - poni_rot1_rad=0, - poni_rot2_rad=0 + poni_rot1_rad=-0.001396263, + poni_rot2_rad=-0.003839724 ) geom = SampleGeometryModel(beam_location_pxl=Coordinate(x=1000,y=1000), @@ -251,6 +251,7 @@ class MainWindow(QMainWindow): self.raster.grid_scan_auto.connect(self.daq.raster_scan_auto) self.data_collection.screening.rotation_scan.connect(self.daq.standard_scan) self.data_collection.simple.rotation_scan.connect(self.daq.standard_scan) + self.data_collection.simple.parameters_changed.connect(self.daq.smart_params) self.raster.grid_scan_size_changed.connect(self.data_collection.raster.grid_scan_size_change) self.status_bar.set_pgroup.connect(self.daq.set_pgroup) diff --git a/gui/src/aaregui/panels/rotation_data_collection.py b/gui/src/aaregui/panels/rotation_data_collection.py index 3e3ac470..ac607235 100644 --- a/gui/src/aaregui/panels/rotation_data_collection.py +++ b/gui/src/aaregui/panels/rotation_data_collection.py @@ -123,7 +123,7 @@ class RotationDataCollectionPanel(ScanSettingsPanel): @Slot() def run_measurement(self): r = RotationScanRequest( - file_prefix=str(add_data_to_path(self._filename)), #self._filename, + file_prefix=str(add_data_to_path(self._filename)), start_omega_deg=self.start_angle.value, steps=self.image_number(), incr_omega_deg=self.image_angle.value, diff --git a/gui/src/aaregui/panels/smart_rotation_panel.py b/gui/src/aaregui/panels/smart_rotation_panel.py index 5c765d87..00d7e3dc 100644 --- a/gui/src/aaregui/panels/smart_rotation_panel.py +++ b/gui/src/aaregui/panels/smart_rotation_panel.py @@ -3,19 +3,22 @@ import math from PySide6.QtCore import Slot, Qt, Signal from PySide6.QtWidgets import QWidget, QGridLayout, QLabel, QFrame, QPushButton, QSpacerItem, QSizePolicy -from aaredaqlib.models import DAQStatusModel +from aaredaqlib.models import DAQStatusModel, SimpleScanParameters from aaredaqlib.rotation_scan import RotationScanRequest +from aaregui.panels.rotation_data_collection import add_data_to_path from aaregui.widgets.number_line_edit import NumberLineEdit class SimpleRotationSettingsPanel(QWidget): rotation_scan = Signal(RotationScanRequest) viewer_track_online = Signal() + parameters_changed = Signal(SimpleScanParameters) def __init__(self, parent=None): super().__init__(parent) + self.n_images = 1 self.xtal_size_dose_rate_MGy_s = None self.xtal_x = None self.xtal_y = None @@ -24,14 +27,12 @@ class SimpleRotationSettingsPanel(QWidget): self.dose_rate_MGy_s = None self._filename = "" self.dtz = 120 - self.n_images = None - self.transmission = 1.0 - self.image_time_s = None - self.image_angle = 0.2 + self.image_time_s = 0.1 self.__d = None self._temperature = 100 self.__omega = 0 - self._wilson_b = 0 + self._wilson_b = None + self.parameters = SimpleScanParameters() self._layout = QGridLayout(self) @@ -44,73 +45,75 @@ class SimpleRotationSettingsPanel(QWidget): self._layout.addWidget(QLabel("Å", parent=self), 0, 4) self.visible_res_enter.newValue.connect(self.set_visible_resolution) - # Angular range (entry) - self._layout.addWidget(QLabel("Total angle", parent=self), 1, 0) - self.angular_range_enter = NumberLineEdit( - 5.0, 1000.0, decimals=3, default=210.0, parent=self - ) - self._layout.addWidget(self.angular_range_enter, 1, 1, 1, 3) + self._layout.addWidget(QLabel("Start angle", parent=self), 1, 0) + self.start_angle_enter = NumberLineEdit(-720, 720.0, 0.0, decimals=3, parent=self) + self._layout.addWidget(self.start_angle_enter, 1, 1, 1, 3) self._layout.addWidget(QLabel("°", parent=self), 1, 4) - self.visible_res_enter.newValue.connect(self.set_visible_resolution) - - # Calculated labels - self._layout.addWidget(QLabel("Target resolution", parent=self), 2, 0) - self.target_res_label = QLabel("--", parent=self) - self.target_res_label.setAlignment(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter) - self._layout.addWidget(self.target_res_label, 2, 1, 1, 3) - self._layout.addWidget(QLabel("Å", parent=self), 2, 4) - - self._layout.addWidget(QLabel("Image angle", parent=self), 3, 0) - self.image_angle_label = QLabel("--", parent=self) - self.image_angle_label.setAlignment(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter) - self._layout.addWidget(self.image_angle_label, 3, 1, 1, 3) - self._layout.addWidget(QLabel("°", parent=self), 3, 4) - - self._layout.addWidget(QLabel("Image time", parent=self), 4, 0) - self.image_time_label = QLabel("--", parent=self) - self.image_time_label.setAlignment(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter) - self._layout.addWidget(self.image_time_label, 4, 1, 1, 3) - self._layout.addWidget(QLabel("s", parent=self), 4, 4) - - self._layout.addWidget(QLabel("Transmission", parent=self), 5, 0) - self.transmission_label = QLabel("--", parent=self) - self.transmission_label.setAlignment(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter) - self._layout.addWidget(self.transmission_label, 5, 1, 1, 3) - self._layout.addWidget(QLabel("%", parent=self), 5, 4) - - self._layout.addWidget(QLabel("Detector distance", parent=self), 6, 0) - self.dtz_label = QLabel(f"--", parent=self) - self.dtz_label.setAlignment(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter) - self._layout.addWidget(self.dtz_label, 6, 1, 1, 3) - self._layout.addWidget(QLabel("mm", parent=self), 6, 4) - - self._layout.addWidget(QLabel("Target Dose", parent=self), 7, 0) - self.target_dose_label = QLabel(f"--", parent=self) - self.target_dose_label.setAlignment(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter) - self._layout.addWidget(self.target_dose_label, 7, 1, 1, 3) - self._layout.addWidget(QLabel("MGy", parent=self), 7, 4) - - self._layout.addWidget(QLabel("Calculated Dose Rate", parent=self), 8, 0) - self.calculated_dose_rate_label = QLabel(f"--", parent=self) - self.calculated_dose_rate_label.setAlignment(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter) - self._layout.addWidget(self.calculated_dose_rate_label, 8, 1, 1, 3) - self._layout.addWidget(QLabel("MGy S^-1", parent=self), 8, 4) - - self._layout.addWidget(QLabel("Temperature", parent=self), 9, 0) - self.temp_enter = NumberLineEdit(80, 330, decimals=2, default=100.0, parent=self) - self._layout.addWidget(self.temp_enter, 9, 1, 1, 3) - self._layout.addWidget(QLabel("K", parent=self), 9, 4) - self.temp_enter.newValue.connect(self.set_temperature) - - self._layout.addWidget(QLabel("Start angle", parent=self), 10, 0) - self.start_angle = NumberLineEdit(-720, 720.0, 0.0, decimals=3, parent=self) - self._layout.addWidget(self.start_angle, 10, 1, 1, 3) - self._layout.addWidget(QLabel("°", parent=self), 10, 4) self.omega_button = QPushButton("Ω") self.omega_button.setFixedWidth(20) self.omega_button.clicked.connect(self.update_omega_start) - self._layout.addWidget(self.omega_button, 10, 5) + self._layout.addWidget(self.omega_button, 1, 5) + + # Angular range (entry) + self._layout.addWidget(QLabel("Total angle", parent=self), 2, 0) + self.angular_range_enter = NumberLineEdit( + 5.0, 1000.0, decimals=3, default=360.0, parent=self + ) + self._layout.addWidget(self.angular_range_enter, 2, 1, 1, 3) + self._layout.addWidget(QLabel("°", parent=self), 2, 4) + self.visible_res_enter.newValue.connect(self.set_total_angle) + + self._layout.addWidget(QLabel("Image angle", parent=self), 3, 0) + self.image_angle_enter = NumberLineEdit( + 0.001, 1.000, decimals=3, default=0.2, parent=self + ) + self._layout.addWidget(self.image_angle_enter, 3, 1, 1, 3) + self._layout.addWidget(QLabel("°", parent=self), 3, 4) + self.image_angle_enter.newValue.connect(self.set_image_angle) + + self._layout.addWidget(QLabel("Temperature", parent=self), 4, 0) + self.temp_enter = NumberLineEdit(80, 330, decimals=2, default=100.0, parent=self) + self._layout.addWidget(self.temp_enter, 4, 1, 1, 3) + self._layout.addWidget(QLabel("K", parent=self), 4, 4) + self.temp_enter.newValue.connect(self.set_temperature) + + # Calculated labels + self._layout.addWidget(QLabel("Target resolution", parent=self), 5, 0) + self.target_res_label = QLabel("--", parent=self) + self.target_res_label.setAlignment(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter) + self._layout.addWidget(self.target_res_label, 5, 1, 1, 3) + self._layout.addWidget(QLabel("Å", parent=self), 5, 4) + + self._layout.addWidget(QLabel("Image time", parent=self), 6, 0) + self.image_time_label = QLabel("--", parent=self) + self.image_time_label.setAlignment(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter) + self._layout.addWidget(self.image_time_label, 6, 1, 1, 3) + self._layout.addWidget(QLabel("s", parent=self), 6, 4) + + self._layout.addWidget(QLabel("Transmission", parent=self), 7, 0) + self.transmission_label = QLabel("--", parent=self) + self.transmission_label.setAlignment(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter) + self._layout.addWidget(self.transmission_label, 7, 1, 1, 3) + self._layout.addWidget(QLabel("%", parent=self), 7, 4) + + self._layout.addWidget(QLabel("Detector distance", parent=self), 8, 0) + self.dtz_label = QLabel(f"--", parent=self) + self.dtz_label.setAlignment(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter) + self._layout.addWidget(self.dtz_label, 8, 1, 1, 3) + self._layout.addWidget(QLabel("mm", parent=self), 8, 4) + + self._layout.addWidget(QLabel("Target Dose", parent=self), 9, 0) + self.target_dose_label = QLabel(f"--", parent=self) + self.target_dose_label.setAlignment(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter) + self._layout.addWidget(self.target_dose_label, 9, 1, 1, 3) + self._layout.addWidget(QLabel("MGy", parent=self), 9, 4) + + self._layout.addWidget(QLabel("Calculated Dose Rate", parent=self), 10, 0) + self.calculated_dose_rate_label = QLabel(f"--", parent=self) + self.calculated_dose_rate_label.setAlignment(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter) + self._layout.addWidget(self.calculated_dose_rate_label, 10, 1, 1, 3) + self._layout.addWidget(QLabel("MGy S^-1", parent=self), 10, 4) self._layout.addWidget(QLabel("Wilson B Factor", parent=self), 11, 0) self.wilson_b_label = QLabel(f"--", parent=self) @@ -178,7 +181,11 @@ class SimpleRotationSettingsPanel(QWidget): @Slot(DAQStatusModel) def update_daq_status(self, s: DAQStatusModel): self.__d = s + #TODO only update best_res after raster finished otherwise ask to update. or have toggle to overwrite with user value + #TODO only take best res from flat face scan + #TODO identify flat face!!!! best_res = s.last_best_res + self.__omega = s.geom.omega_deg if best_res is not None: # clamp to control limits and update field; this will also trigger recalculation @@ -187,10 +194,10 @@ class SimpleRotationSettingsPanel(QWidget): self.visible_res_enter.update_value(v) self.set_visible_resolution(v) - best_b_factor = s.last_best_b_factor + self._wilson_b = s.last_best_b_factor - if best_b_factor is not None: - self.wilson_b_label.setText(f"{best_b_factor:.2f}") + if self._wilson_b is not None: + self.wilson_b_label.setText(f"{self._wilson_b:.2f}") xtal_size = s.crystal_size if xtal_size is not None: @@ -211,12 +218,21 @@ class SimpleRotationSettingsPanel(QWidget): self.wavelength_label.setText("N/A") else: self.wavelength_label.setText(f"{s.diffraction.wavelength_angstrom:.3f}") + self.update_calculated_labels() @Slot(float) def set_visible_resolution(self, v: float): self.update_calculated_labels() + @Slot(float) + def set_total_angle(self, v: float): + self.update_calculated_labels() + + @Slot(float) + def set_image_angle(self, v: float): + self.update_calculated_labels() + @Slot(float) def set_temperature(self, v: float): self.update_calculated_labels() @@ -227,7 +243,7 @@ class SimpleRotationSettingsPanel(QWidget): @Slot() def update_omega_start(self): - self.start_angle.update_value(self.__omega) + self.start_angle_enter.update_value(self.__omega) def update_calculated_labels(self): if self.__d is None: @@ -238,9 +254,12 @@ class SimpleRotationSettingsPanel(QWidget): #TODO read resolution estiamtion from jfjoch total_angle = self.angular_range_enter.value d_vis = self.visible_res_enter.value + image_angle = self.image_angle_enter.value + if d_vis == 0.0: + d_vis = 1.3 d_tar = 1/(1/d_vis + 0.1) + self.target_res_label.setText(f"{d_tar:.2f}") - #self.wilson_b_label.setText("--") Kdose = 2000 / (self.__d.diffraction.wavelength_angstrom**2) @@ -271,10 +290,9 @@ class SimpleRotationSettingsPanel(QWidget): total_time_s = self.target_dose_MGy / self.dose_rate_MGy_s self.calculated_dose_label.setText(f"{self.xtal_size_dose_rate_MGy_s*total_time_s:.2f}") - - self.image_angle = 0.2 - self.image_angle_label.setText(f"{self.image_angle:.3f}") - self.n_images = total_angle / self.image_angle + if image_angle == 0.0: + image_angle = 0.001 + self.n_images = round(total_angle / image_angle) self.image_time_s = total_time_s / self.n_images if self.image_time_s < 0.01: self.transmission = self.image_time_s / 0.01 @@ -295,13 +313,23 @@ class SimpleRotationSettingsPanel(QWidget): else: self.dtz_label.setText(f"{self.dtz:.2f}") + self.parameters = SimpleScanParameters( + dtz=int(round(self.dtz)), + exp_time_s=self.image_time_s, + start_omega_deg=self.start_angle_enter.value, + incr_omega_deg=image_angle, + steps=self.n_images, + transmission=self.transmission + ) + self.parameters_changed.emit(self.parameters) + @Slot() def run_measurement(self): r = RotationScanRequest( - file_prefix=self._filename, - start_omega_deg=0,#,self.start_angle.value, + file_prefix=str(add_data_to_path(self._filename)), + start_omega_deg=self.start_angle_enter.value, steps=self.n_images, - incr_omega_deg=self.image_angle, + incr_omega_deg=self.image_angle_enter.value, dtz=self.dtz, transmission=self.transmission, screening=False, diff --git a/gui/src/aaregui/threads/daq_worker.py b/gui/src/aaregui/threads/daq_worker.py index 3656cdde..ee69b9b6 100644 --- a/gui/src/aaregui/threads/daq_worker.py +++ b/gui/src/aaregui/threads/daq_worker.py @@ -7,7 +7,7 @@ from jfjoch_client import ScanResult, ScanResultImagesInner from aaredaqlib.coordinate import SmargonCoordinate, Coordinate from aaredaqlib.models import DAQStatusModel, SampleShortInfoList, SampleShortInfo, SampleCameraSettings, \ - AutofocusSettings + AutofocusSettings, SimpleScanParameters from aaredaqlib.raster_grid import RasterGridRequest, CompletedRasterGrid from aaredaqlib.rotation_scan import RotationScanRequest, CompletedRotationScan @@ -351,6 +351,13 @@ class DAQWorker(QObject): reply = self.__net_manager.post(request, QByteArray(body.encode("utf-8"))) reply.finished.connect(lambda: self.handle_auto_scan_response(reply, s.db_id)) + @Slot(SimpleScanParameters) + def smart_params(self, p: SimpleScanParameters): + if self.__base_url is None: + print(f"POST /scan/smart_params: {p.model_dump_json()}") + return + self.generic_post("scan/smart_params", p.model_dump_json()) + @Slot(Coordinate) def abr_tweak(self, c: Coordinate): self.generic_post("beamline/tweak_abr_meas_pos", c.model_dump_json())