From be3fd302296083690bebd1bad7cc30b58359d468 Mon Sep 17 00:00:00 2001 From: ahmedk <karim.ahmed@xfel.eu> Date: Wed, 30 Jun 2021 15:59:13 +0200 Subject: [PATCH] split from epix changes and refactor --- ...Jungfrau_Gain_Correct_and_Verify_NBC.ipynb | 4 +- ...rk_analysis_all_gains_burst_mode_NBC.ipynb | 62 ++--- .../Characterize_Darks_ePix100_NBC.ipynb | 178 ++++++------ .../ePix100/Correction_ePix100_NBC.ipynb | 263 ++++++++++-------- 4 files changed, 254 insertions(+), 253 deletions(-) diff --git a/notebooks/Jungfrau/Jungfrau_Gain_Correct_and_Verify_NBC.ipynb b/notebooks/Jungfrau/Jungfrau_Gain_Correct_and_Verify_NBC.ipynb index a9b0c59a9..3e5cbd264 100644 --- a/notebooks/Jungfrau/Jungfrau_Gain_Correct_and_Verify_NBC.ipynb +++ b/notebooks/Jungfrau/Jungfrau_Gain_Correct_and_Verify_NBC.ipynb @@ -45,7 +45,7 @@ "chunk_size_idim = 1 # chunking size of imaging dimension, adjust if user software is sensitive to this.\n", "integration_time = 4.96 # integration time in us, will be overwritten by value in file\n", "mem_cells = 0 # leave memory cells equal 0, as it is saved in control information starting 2019.\n", - "db_module = [\"\"] # ID of module in calibration database\n", + "db_module = \"\" # ID of module in calibration database\n", "manual_slow_data = False # if true, use manually entered bias_voltage and integration_time values\n", "chunk_size = 0\n", "\n", @@ -66,7 +66,6 @@ "import time\n", "import traceback\n", "import warnings\n", - "from extra_data import RunDirectory\n", "from functools import partial\n", "from pathlib import Path\n", "\n", @@ -76,6 +75,7 @@ "import numpy as np\n", "import tabulate\n", "from IPython.display import Latex, Markdown, display\n", + "from extra_data import RunDirectory\n", "from matplotlib.colors import LogNorm\n", "\n", "from cal_tools import jungfraulib\n", diff --git a/notebooks/Jungfrau/Jungfrau_dark_analysis_all_gains_burst_mode_NBC.ipynb b/notebooks/Jungfrau/Jungfrau_dark_analysis_all_gains_burst_mode_NBC.ipynb index d628865ec..337a7a67c 100644 --- a/notebooks/Jungfrau/Jungfrau_dark_analysis_all_gains_burst_mode_NBC.ipynb +++ b/notebooks/Jungfrau/Jungfrau_dark_analysis_all_gains_burst_mode_NBC.ipynb @@ -139,38 +139,6 @@ "report = get_report(out_folder)" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def check_memory_cells(\n", - " run_dir,\n", - " karabo_id_control,\n", - " manual_slow_data,\n", - " integration_time,\n", - " bias_voltage,\n", - "):\n", - "\n", - " n_scs = int(run_dir.get_array(\n", - " f\"{karabo_id_control}/DET/CONTROL\", \"storageCells.value\")[0]) + 1\n", - " sc_s = int(run_dir.get_array(\n", - " f\"{karabo_id_control}/DET/CONTROL\", \"storageCellStart.value\")[0])\n", - " # TODO: Read RUN source as previously done?\n", - " # TODO: Is this still needed?\n", - " if not manual_slow_data:\n", - " integration_time = float(run_dir.get_array(\n", - " f\"{karabo_id_control}/DET/CONTROL\", \"exposureTime.value\")[0]) * 1e6\n", - " bias_voltage = int(np.squeeze(run_dir.get_array(\n", - " f\"{karabo_id_control}/DET/CONTROL\", \"vHighVoltage.value\"))[0, 0])\n", - "\n", - " print(f\"Integration time is {integration_time} us\")\n", - " print(f\"Bias voltage is {bias_voltage} V\")\n", - "\n", - " return n_scs, sc_s, integration_time, bias_voltage" - ] - }, { "cell_type": "code", "execution_count": null, @@ -180,6 +148,8 @@ "noise_map = OrderedDict()\n", "offset_map = OrderedDict()\n", "memory_cells = None\n", + "\n", + "\n", "for mod in karabo_da:\n", "\n", " h5_inst_path = h5path.format(\n", @@ -197,23 +167,27 @@ " n_empty_sc = 0\n", "\n", " run_dir = RunDirectory(f\"{in_folder}/r{r_n:04d}/\")\n", - "\n", - " run_mcells, sc_start, integration_time, bias_voltage = jungfraulib.check_memory_cells( # noqa\n", - " run_dir,\n", - " karabo_id_control,\n", - " manual_slow_data,\n", - " integration_time,\n", - " bias_voltage,\n", - " )\n", - "\n", + " \n", " if mod not in noise_map.keys():\n", "\n", + " (run_mcells, sc_start, integration_time, bias_voltage) = jungfraulib.check_memory_cells( # noqa\n", + " run_dir,\n", + " karabo_id_control,\n", + " manual_slow_data,\n", + " integration_time,\n", + " bias_voltage,\n", + " )\n", + " print(f\"Integration time is {integration_time} us\")\n", + " print(f\"Bias voltage is {bias_voltage} V\")\n", + "\n", " if run_mcells == 1:\n", " memory_cells = 1\n", - " print('Dark runs in single cell mode\\n storage cell start: {:02d}'.format(sc_start))\n", + " print('Dark runs in single cell mode\\n'\n", + " f'\\nstorage cell start: {sc_start:02d}')\n", " else:\n", " memory_cells = 16\n", - " print('Dark runs in burst mode\\n storage cell start: {:02d}'.format(sc_start))\n", + " print('Dark runs in burst mode\\n'\n", + " f'storage cell start: {sc_start:02d}')\n", "\n", " noise_map[mod] = np.zeros(sensor_size+[memory_cells, 3])\n", " offset_map[mod] = np.zeros(sensor_size+[memory_cells, 3])\n", @@ -403,7 +377,7 @@ " mdn = np.nanmedian(d, axis=(0, 1))[None, None, :, :]\n", " std = np.nanstd(d, axis=(0, 1))[None, None, :, :] \n", " idx = (d > badpixel_threshold_sigma*std+mdn) | (d < (-badpixel_threshold_sigma)*std+mdn)\n", - " \n", + "\n", " return idx" ] }, diff --git a/notebooks/ePix100/Characterize_Darks_ePix100_NBC.ipynb b/notebooks/ePix100/Characterize_Darks_ePix100_NBC.ipynb index 5ccf0f35c..89ea1f460 100644 --- a/notebooks/ePix100/Characterize_Darks_ePix100_NBC.ipynb +++ b/notebooks/ePix100/Characterize_Darks_ePix100_NBC.ipynb @@ -30,7 +30,7 @@ "receiver_id = \"RECEIVER\" # inset for receiver devices\n", "path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # the template to use to access data\n", "h5path = '/INSTRUMENT/{}/DET/{}:daqOutput/data/image/pixels' # path in the HDF5 file to images\n", - "h5path_t = 'data.backTemp' # path to find temperature at\n", + "h5path_t = '/INSTRUMENT/{}/DET/{}:daqOutput/data/backTemp' # path to find temperature at\n", "h5path_cntrl = '/CONTROL/{}/DET' # path to control data\n", "\n", "use_dir_creation_date = True\n", @@ -45,9 +45,7 @@ "bias_voltage = 200 # bias voltage\n", "in_vacuum = False # detector operated in vacuum\n", "fix_temperature = 290. # fix temperature to this value\n", - "operation_mode = '' # Detector operation mode, optional\n", - "min_trains = 1 # Minimum number of trains that should be available to process dark constants. Default 1.\n", - "max_trains = 500 # Maximum number of trains to use for processing dark constants. Set to 0 to use all available trains." + "operation_mode = '' # Detector operation mode, optional" ] }, { @@ -59,14 +57,14 @@ "import os\n", "import warnings\n", "\n", - "import numpy as np\n", - "from extra_data import RunDirectory\n", + "warnings.filterwarnings('ignore')\n", "\n", - "import XFELDetAna.xfelprofiler as xprof\n", - "from XFELDetAna import xfelpyanatools as xana\n", - "from XFELDetAna import xfelpycaltools as xcal\n", - "from XFELDetAna.plotting.util import prettyPlotting\n", - "from XFELDetAna.util import env\n", + "import h5py\n", + "import matplotlib.pyplot as plt\n", + "from IPython.display import Latex, Markdown, display\n", + "\n", + "%matplotlib inline\n", + "import numpy as np\n", "from cal_tools.tools import (\n", " get_dir_creation_date,\n", " get_pdu_from_db,\n", @@ -75,27 +73,32 @@ " save_const_to_h5,\n", " send_to_db,\n", ")\n", - "from iCalibrationDB import Conditions, Constants" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "env.iprofile = cluster_profile\n", + "from iCalibrationDB import Conditions, Constants, Detectors, Versions\n", + "from iCalibrationDB.detectors import DetectorTypes\n", + "from XFELDetAna.util import env\n", "\n", - "warnings.filterwarnings('ignore')\n", + "env.iprofile = cluster_profile\n", + "from XFELDetAna import xfelpyanatools as xana\n", + "from XFELDetAna import xfelpycaltools as xcal\n", + "from XFELDetAna.detectors.fastccd import readerh5 as fastccdreaderh5\n", + "from XFELDetAna.plotting.util import prettyPlotting\n", "\n", "prettyPlotting = True\n", + "import XFELDetAna.xfelprofiler as xprof\n", "\n", "profiler = xprof.Profiler()\n", "profiler.disable()\n", + "from XFELDetAna.xfelreaders import ChunkReader\n", "\n", "h5path = h5path.format(karabo_id, receiver_id)\n", "h5path_t = h5path_t.format(karabo_id, receiver_id)\n", - "h5path_cntrl = h5path_cntrl.format(karabo_id)" + "h5path_cntrl = h5path_cntrl.format(karabo_id)\n", + "\n", + "def nImagesOrLimit(nImages, limit):\n", + " if limit == 0:\n", + " return nImages\n", + " else:\n", + " return min(nImages, limit)" ] }, { @@ -110,15 +113,14 @@ "\n", "report = get_report(out_folder)\n", "\n", - "x = 708 # rows of the ePix100\n", - "y = 768 # columns of the ePix100\n", - "run_parallel = False\n", + "x = 708 # rows of the xPix100\n", + "y = 768 # columns of the xPix100\n", "\n", "ped_dir = os.path.join(in_folder, f\"r{run:04d}\")\n", "fp_name = path_template.format(run, karabo_da[0]).format(sequence)\n", - "run_dir = RunDirectory(ped_dir)\n", + "filename = os.path.join(ped_dir, fp_name)\n", "\n", - "print(f\"Reading data from: {os.path.join(ped_dir, fp_name)}\\n\")\n", + "print(f\"Reading data from: {filename}\\n\")\n", "print(f\"Run number: {run}\")\n", "print(f\"HDF5 path: {h5path}\")\n", "if use_dir_creation_date:\n", @@ -134,49 +136,31 @@ "outputs": [], "source": [ "sensorSize = [x, y]\n", - "chunk_size = 100 # Number of images to read per chunk.\n", + "chunkSize = 100 #Number of images to read per chunk\n", "\n", - "# Sensor area will be analysed according to blocksize\n", + "#Sensor area will be analysed according to blocksize\n", "blockSize = [sensorSize[0] // 2, sensorSize[1] // 2]\n", "xcal.defaultBlockSize = blockSize\n", "cpuCores = 4 #Specifies the number of running cpu cores\n", "memoryCells = 1 #No mamery cells\n", "\n", - "h5_data = f\"{karabo_id}/DET/{receiver_id}:daqOutput\"\n", - "\n", - "# Specifies total number of images to proceed\n", - "n_imgs = run_dir.get_data_counts(\n", - " h5_data, \"data.image.pixels\").shape[0]\n", - "\n", - "# Modify n_imgs to process based on given maximum\n", - "# and minimun number of trains.\n", - "if max_trains and n_imgs > max_trains:\n", - " n_imgs = max_trains\n", - "\n", - "if n_imgs < min_trains:\n", - " raise ValueError(\n", - " f\"Less than {min_trains} trains are available in RAW data.\"\n", - " \" Not enough data to process darks.\")\n", - "\n", - "print(f\"Number of dark images to analyze: {n_imgs}.\")\n", - "\n", - "integration_time = int(run_dir.get_array(\n", - " f\"{karabo_id}/DET/CONTROL\",\n", - " \"expTime.value\")[0])\n", - "temperature = np.mean(run_dir.get_array(\n", - " f\"{karabo_id}/DET/{receiver_id}:daqOutput\",\n", - " f\"{h5path_t}\").values) / 100.\n", + "#Specifies total number of images to proceed\n", + "nImages = fastccdreaderh5.getDataSize(filename, h5path)[0]\n", + "nImages = nImagesOrLimit(nImages, number_dark_frames)\n", + "print(\"\\nNumber of dark images to analyze: \", nImages)\n", + "run_parallel = False\n", "\n", - "if fix_temperature != 0:\n", - " temperature_k = fix_temperature\n", - " print(\"Temperature is fixed!\")\n", - "else:\n", + "with h5py.File(filename, 'r') as f:\n", + " integration_time = int(f[os.path.join(h5path_cntrl, \"CONTROL\",\"expTime\", \"value\")][0])\n", + " temperature = np.mean(f[h5path_t])/100.\n", " temperature_k = temperature + 273.15\n", - "\n", - "print(f\"Bias voltage is {bias_voltage} V\")\n", - "print(f\"Detector integration time is set to {integration_time}\")\n", - "print(f\"Mean temperature was {temperature:0.2f} °C / {temperature_k:0.2f} K\")\n", - "print(f\"Operated in vacuum: {in_vacuum} \")" + " if fix_temperature != 0:\n", + " temperature_k = fix_temperature\n", + " print(\"Temperature is fixed!\")\n", + " print(f\"Bias voltage is {bias_voltage} V\")\n", + " print(f\"Detector integration time is set to {integration_time}\")\n", + " print(f\"Mean temperature was {temperature:0.2f} °C / {temperature_k:0.2f} K\")\n", + " print(f\"Operated in vacuum: {in_vacuum} \")" ] }, { @@ -185,22 +169,11 @@ "metadata": {}, "outputs": [], "source": [ - "noiseCal = xcal.NoiseCalculator(\n", - " sensorSize, memoryCells,\n", - " cores=cpuCores,\n", - " blockSize=blockSize,\n", - " parallel=run_parallel,\n", - ")\n", - "\n", - "for chunk in run_dir.select_trains(np.s_[:n_imgs]).split_trains(n_imgs//chunk_size):\n", - "\n", - " data = chunk.get_array(h5_data, \"data.image.pixels\")\n", - " data = np.moveaxis(data.values.astype(np.float), 0, 2)\n", - " dx = np.count_nonzero(data, axis=(0, 1))\n", - " data = data[:, :, dx != 0]\n", - "\n", - " # Fill calculators with data\n", - " noiseCal.fill(data)" + "reader = ChunkReader(filename, fastccdreaderh5.readData,\n", + " nImages, chunkSize,\n", + " path=h5path,\n", + " pixels_x=sensorSize[0],\n", + " pixels_y=sensorSize[1], )" ] }, { @@ -209,12 +182,33 @@ "metadata": {}, "outputs": [], "source": [ + "noiseCal = xcal.NoiseCalculator(sensorSize, memoryCells,\n", + " cores=cpuCores, blockSize=blockSize,\n", + " parallel=run_parallel)\n", + "histCalRaw = xcal.HistogramCalculator(sensorSize, bins=1000,\n", + " range=[0, 10000], parallel=False,\n", + " memoryCells=memoryCells,\n", + " cores=cpuCores, blockSize=blockSize)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for data in reader.readChunks():\n", + " dx = np.count_nonzero(data, axis=(0, 1))\n", + " data = data[:, :, dx != 0]\n", + " histCalRaw.fill(data)\n", + " noiseCal.fill(data) #Fill calculators with data\n", + " \n", "constant_maps = {}\n", - "constant_maps['Offset'] = noiseCal.getOffset() # Produce offset map\n", - "constant_maps['Noise'] = noiseCal.get() # Produce noise map\n", + "constant_maps['Offset'] = noiseCal.getOffset() #Produce offset map\n", + "constant_maps['Noise'] = noiseCal.get() #Produce noise map\n", "\n", - "noiseCal.reset() # Reset noise calculator\n", - "print(\"Initial constant maps are created\")" + "noiseCal.reset() #Reset noise calculator\n", + "print(\"Initial maps were created\")" ] }, { @@ -270,13 +264,18 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": false + }, "outputs": [], "source": [ "# Save constants to DB\n", "dclass=\"ePix100\"\n", "md = None\n", "\n", + "# If PDU(db_module) is not given with input parameters \n", + "# retrieve the connected physical detector unit\n", + "\n", "for const_name in constant_maps.keys():\n", " det = getattr(Constants, dclass)\n", " const = getattr(det, const_name)()\n", @@ -292,6 +291,10 @@ " parm.lower_deviation = temp_limits\n", " parm.upper_deviation = temp_limits\n", "\n", + " # This should be used in case of running notebook \n", + " # by a different method other than myMDC which already\n", + " # sends CalCat info.\n", + " # TODO: Set db_module to \"\" by default in the first cell\n", " if not db_module:\n", " db_module = get_pdu_from_db(karabo_id, karabo_da, const,\n", " condition, cal_db_interface,\n", @@ -314,6 +317,13 @@ " f\"• Temperature: {temperature_k}\\n• In Vacuum: {in_vacuum}\\n\"\n", " f\"• Creation time: {md.calibration_constant_version.begin_at if md is not None else creation_time}\\n\")" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/notebooks/ePix100/Correction_ePix100_NBC.ipynb b/notebooks/ePix100/Correction_ePix100_NBC.ipynb index 8c4e9a9a2..0c46331ff 100644 --- a/notebooks/ePix100/Correction_ePix100_NBC.ipynb +++ b/notebooks/ePix100/Correction_ePix100_NBC.ipynb @@ -25,12 +25,12 @@ "\n", "karabo_id = \"MID_EXP_EPIX-1\" # karabo karabo_id\n", "karabo_da = \"EPIX01\" # data aggregators\n", - "db_module = \"\" # module id in the database\n", + "db_module = \"ePix100_M15\" # module id in the database\n", "receiver_id = \"RECEIVER\" # inset for receiver devices\n", "path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # the template to use to access data\n", - "h5path = '{}/DET/{}:daqOutput' # path in the HDF5 file to images\n", - "h5path_t = 'data.backTemp' # path to find temperature at # TODO: modify in calibration_configurations.\n", - "h5path_cntrl = '/CONTROL/{}/DET' # path to control data # TODO: remove from calibration_configurations.\n", + "h5path = '/INSTRUMENT/{}/DET/{}:daqOutput/data/image' # path in the HDF5 file to images\n", + "h5path_t = '/INSTRUMENT/{}/DET/{}:daqOutput/data/backTemp' # path to find temperature at\n", + "h5path_cntrl = '/CONTROL/{}/DET' # path to control data\n", "\n", "use_dir_creation_date = True # date constants injected before directory creation time\n", "cal_db_interface = \"tcp://max-exfl016:8015#8025\" # calibration DB interface to use\n", @@ -50,9 +50,8 @@ "\n", "relative_gain = False # Apply relative gain correction.\n", "common_mode = True # Apply common mode correction.\n", - "cm_min_frac = 0.25 # No CM correction is performed if after masking the ratio of good pixels falls below this \n", - "cm_noise_sigma = 5. # CM correction noise standard deviation\n", - "temp_deviations = 5. # temperature deviation for the constant operating conditions\n", + "cm_min_frac = 0.25 # No CM correction is performed if after masking the ratio of good pixels falls below this \n", + "cm_noise_sigma = 5. # CM correction noise standard deviation\n", "\n", "split_evt_primary_threshold = 7. # primary threshold for split event correction\n", "split_evt_secondary_threshold = 5. # secondary threshold for split event correction\n", @@ -75,9 +74,7 @@ "\n", "import h5py\n", "import numpy as np\n", - "from IPython.display import Latex, Markdown, display\n", - "\n", - "from extra_data import RunDirectory\n", + "from IPython.display import Latex, display\n", "from pathlib import Path\n", "\n", "import XFELDetAna.xfelprofiler as xprof\n", @@ -115,6 +112,8 @@ "pattern_classification = False # do clustering.\n", "\n", "h5path = h5path.format(karabo_id, receiver_id)\n", + "h5path_t = h5path_t.format(karabo_id, receiver_id)\n", + "h5path_cntrl = h5path_cntrl.format(karabo_id)\n", "plot_unit = 'ADU'\n", "\n", "if relative_gain:\n", @@ -134,7 +133,6 @@ "\n", "in_folder = Path(in_folder)\n", "ped_dir = in_folder / f\"r{run:04d}\"\n", - "run_dir = RunDirectory(ped_dir)\n", "fp_name = path_template.format(run, karabo_da)\n", "\n", "print(f\"Reading from: {ped_dir / fp_name}\")\n", @@ -166,27 +164,21 @@ "# Read slow data from the first available sequence file.\n", "filename = ped_dir / fp_name.format(\n", " sequences[0] if sequences[0] != -1 else 0)\n", - "\n", - "integration_time = int(run_dir.get_array(\n", - " f\"{karabo_id}/DET/CONTROL\",\n", - " \"expTime.value\")[0])\n", - "temperature = np.mean(run_dir.get_array(\n", - " f\"{karabo_id}/DET/{receiver_id}:daqOutput\",\n", - " f\"{h5path_t}\").values) / 100.\n", - "\n", - "if fix_temperature != 0:\n", - " temperature_k = fix_temperature\n", - " print(\"Temperature is fixed!\")\n", - "else:\n", + "with h5py.File(filename, 'r') as f:\n", + " integration_time = int(\n", + " f[f\"{h5path_cntrl}/CONTROL/expTime/value\"][0])\n", + " temperature = np.mean(f[h5path_t]) / 100.\n", " temperature_k = temperature + 273.15\n", - "\n", - "print(f\"Bias voltage is {bias_voltage} V\")\n", - "print(f\"Detector integration time is set to {integration_time}\")\n", - "print(\n", - " f\"Mean temperature was {temperature:0.2f} °C \"\n", - " f\"/ {temperature_k:0.2f} K at beginning of the run.\"\n", - ")\n", - "print(f\"Operated in vacuum: {in_vacuum} \")\n", + " if fix_temperature != 0:\n", + " temperature_k = fix_temperature\n", + " print(\"Temperature is fixed!\")\n", + " print(f\"Bias voltage is {bias_voltage} V\")\n", + " print(f\"Detector integration time is set to {integration_time}\")\n", + " print(\n", + " f\"Mean temperature was {temperature:0.2f} °C \"\n", + " f\"/ {temperature_k:0.2f} K at beginning of run\"\n", + " )\n", + " print(f\"Operated in vacuum: {in_vacuum} \")\n", "\n", "out_folder = Path(out_folder)\n", "if out_folder.is_dir() and not overwrite:\n", @@ -243,6 +235,8 @@ "metadata": {}, "outputs": [], "source": [ + "temp_limits = 5.\n", + "\n", "cond_dict = {\n", " \"bias_voltage\": bias_voltage,\n", " \"integration_time\": integration_time,\n", @@ -273,8 +267,8 @@ " # TODO: Fix this logic.\n", " for parm in condition.parameters:\n", " if parm.name == \"Sensor Temperature\":\n", - " parm.lower_deviation = temp_deviations\n", - " parm.upper_deviation = temp_deviations\n", + " parm.lower_deviation = temp_limits\n", + " parm.upper_deviation = temp_limits\n", "\n", " const_data[cname] = get_constant_from_db(\n", " karabo_id=karabo_id,\n", @@ -289,11 +283,10 @@ " )\n", "\n", "if relative_gain and const_data[\"RelativeGain\"] is None:\n", - " display(Markdown(\n", - " \"<span style='color: #ffcc00'>WARNING: </span>\" \n", + " print(\n", " \"WARNING: RelativeGain map is requested, but not found./n\"\n", - " \"No gain correction will be applied.\"\n", - " ))\n", + " \"No gain correction will be applied\"\n", + " )\n", " relative_gain = False\n", " plot_unit = 'ADU'" ] @@ -474,94 +467,118 @@ "metadata": {}, "outputs": [], "source": [ - "# Copy and sanitize raw data into corrected data,\n", - "# without the data.image.pixels.\n", - "out_file = out_folder / seq_files[0].name.replace(\"RAW\", \"CORR\")\n", - "\n", - "run_dir.select(\n", - " f\"*{karabo_id}*\", \"*\").deselect(h5path, \"data.image.pixels\").write(out_file) # noqa\n", - "\n", - "with h5py.File(out_file, \"r+\") as ofile:\n", - " try:\n", - " data = run_dir.get_array(\n", - " h5path, \"data.image.pixels\").values\n", - " data = np.compress(\n", - " np.any(data > 0, axis=(1, 2)), data, axis=0)\n", - "\n", - " if limit_images > 0:\n", - " data = data[:limit_images, ...]\n", - "\n", - " oshape = data.shape\n", - " data = np.moveaxis(data, 0, 2)\n", - " # Create dataset in CORR h5 file and add corrected images.\n", - " ddset = ofile.create_dataset(\n", - " \"/INSTRUMENT/\"+h5path+\"/data/image/pixels\",\n", - " oshape,\n", - " chunks=(chunk_size_idim, oshape[1], oshape[2]),\n", - " dtype=np.float32)\n", - "\n", - " # Offset correction.\n", - " data = offsetCorrection.correct(data.astype(np.float32))\n", - " # Common Mode correction.\n", - " if common_mode:\n", - " # Block CM\n", - " data = cmCorrectionB.correct(data)\n", - " # Row CM\n", - " data = cmCorrectionR.correct(data)\n", - " # COL CM\n", - " data = cmCorrectionC.correct(data)\n", - "\n", - " histCalCMCor.fill(data)\n", - " # relative gain correction.\n", - " if relative_gain:\n", - " data = gainCorrection.correct(data.astype(np.float32))\n", - " if photon_energy > 0:\n", - " data /= photon_energy\n", - "\n", - " histCalOffsetCor.fill(data)\n", - " ddset[...] = np.moveaxis(data, 2, 0)\n", - "\n", - " \"\"\"The gain correction is currently applying\n", - " an absolute correction (not a relative correction\n", - " as the implied by the name);\n", - " it changes the scale (the unit of measurement)\n", - " of the data from ADU to either keV or n_of_photons.\n", - " But the pattern classification relies on comparing\n", - " data with the noise map, which is still in ADU.\n", - "\n", - " The best solution is to do a relative gain\n", - " correction first and apply the global absolute\n", - " gain to the data at the end, after clustering.\n", - " \"\"\"\n", - "\n", - " # TODO: Fix conflict between pattern classification\n", - " # and gain corr.\n", - " if pattern_classification:\n", - " ddsetc = ofile.create_dataset(\n", - " h5path+\"/pixels_classified\",\n", - " oshape,\n", - " chunks=(chunk_size_idim, oshape[1], oshape[2]),\n", - " dtype=np.float32, compression=\"gzip\")\n", - "\n", - " ddsetp = ofile.create_dataset(\n", - " h5path+\"/patterns\",\n", + "def copy_and_sanitize_non_cal_data(\n", + " infile: h5py,\n", + " outfile: h5py,\n", + " h5base: str\n", + "):\n", + " \"\"\" Copy and sanitize data in `infile`\n", + " that is not touched by `correctEPIX`. \"\"\"\n", + "\n", + " if h5base.startswith(\"/\"):\n", + " h5base = h5base[1:]\n", + " dont_copy = [h5base+\"/pixels\"]\n", + "\n", + " def visitor(k, item):\n", + " if k not in dont_copy:\n", + " if isinstance(item, h5py.Group):\n", + " outfile.create_group(k)\n", + " elif isinstance(item, h5py.Dataset):\n", + " group = str(k).split(\"/\")\n", + " group = \"/\".join(group[:-1])\n", + " infile.copy(k, outfile[group])\n", + "\n", + " infile.visititems(visitor)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for f in seq_files:\n", + " data = None\n", + " out_file = out_folder / f.name.replace(\"RAW\", \"CORR\")\n", + " with h5py.File(f, \"r\") as infile, h5py.File(out_file, \"w\") as ofile: # noqa\n", + " try:\n", + " copy_and_sanitize_non_cal_data(infile, ofile, h5path)\n", + " data = infile[h5path+\"/pixels\"][()]\n", + " data = np.compress(\n", + " np.any(data > 0, axis=(1, 2)), data, axis=0)\n", + " if limit_images > 0:\n", + " data = data[:limit_images, ...]\n", + "\n", + " oshape = data.shape\n", + " data = np.moveaxis(data, 0, 2)\n", + " ddset = ofile.create_dataset(\n", + " h5path+\"/pixels\",\n", " oshape,\n", " chunks=(chunk_size_idim, oshape[1], oshape[2]),\n", - " dtype=np.int32, compression=\"gzip\")\n", - "\n", - "\n", - " data, patterns = patternClassifier.classify(data)\n", - "\n", - " data[data < (split_evt_primary_threshold*const_data[\"Noise\"])] = 0 # noqa\n", - " ddsetc[...] = np.moveaxis(data, 2, 0)\n", - " ddsetp[...] = np.moveaxis(patterns, 2, 0)\n", - "\n", - " data[patterns != 100] = np.nan\n", - " histCalSECor.fill(data)\n", - " except Exception as e:\n", - " display(Markdown(\n", - " \"<span style='color: #ff0000'>ERROR \"\n", - " f\"</span> While correcting data : {e} .\"))" + " dtype=np.float32)\n", + "\n", + " # Offset correction.\n", + " data = offsetCorrection.correct(data.astype(np.float32))\n", + "\n", + " # Common Mode correction.\n", + " if common_mode:\n", + " # Block CM\n", + " data = cmCorrectionB.correct(data)\n", + " # Row CM\n", + " data = cmCorrectionR.correct(data)\n", + " # COL CM\n", + " data = cmCorrectionC.correct(data)\n", + "\n", + " histCalCMCor.fill(data)\n", + "\n", + " # relative gain correction.\n", + " if relative_gain:\n", + " data = gainCorrection.correct(data.astype(np.float32))\n", + " if photon_energy > 0:\n", + " data /= photon_energy\n", + "\n", + " histCalOffsetCor.fill(data)\n", + " ddset[...] = np.moveaxis(data, 2, 0)\n", + "\n", + " \"\"\"The gain correction is currently applying\n", + " an absolute correction (not a relative correction\n", + " as the implied by the name);\n", + " it changes the scale (the unit of measurement)\n", + " of the data from ADU to either keV or n_of_photons.\n", + " But the pattern classification relies on comparing\n", + " data with the noise map, which is still in ADU.\n", + "\n", + " The best solution is to do a relative gain\n", + " correction first and apply the global absolute\n", + " gain to the data at the end, after clustering.\n", + " \"\"\"\n", + " \n", + " # TODO: Fix conflict between pattern classification\n", + " # and gain corr.\n", + " if pattern_classification:\n", + " ddsetc = ofile.create_dataset(\n", + " h5path+\"/pixels_classified\",\n", + " oshape,\n", + " chunks=(chunk_size_idim, oshape[1], oshape[2]),\n", + " dtype=np.float32, compression=\"gzip\")\n", + "\n", + " ddsetp = ofile.create_dataset(\n", + " h5path+\"/patterns\",\n", + " oshape,\n", + " chunks=(chunk_size_idim, oshape[1], oshape[2]),\n", + " dtype=np.int32, compression=\"gzip\")\n", + "\n", + "\n", + " data, patterns = patternClassifier.classify(data)\n", + "\n", + " data[data < (split_evt_primary_threshold*const_data[\"Noise\"])] = 0 # noqa\n", + " ddsetc[...] = np.moveaxis(data, 2, 0)\n", + " ddsetp[...] = np.moveaxis(patterns, 2, 0)\n", + "\n", + " data[patterns != 100] = np.nan\n", + " histCalSECor.fill(data)\n", + " except Exception as e:\n", + " print(f\"ERROR applying common mode correction for {f}: {e}\")" ] }, { -- GitLab