diff --git a/notebooks/Jungfrau/create_gain_map.ipynb b/notebooks/Jungfrau/create_gain_map.ipynb
index 4979af6d12af715c47e33c3a35d467618ad76270..5583c6608f16c313a0b52427dd018389fbceb890 100755
--- a/notebooks/Jungfrau/create_gain_map.ipynb
+++ b/notebooks/Jungfrau/create_gain_map.ipynb
@@ -25,14 +25,14 @@
"outputs": [],
"source": [
"in_folder = '/gpfs/exfel/exp/MID/202330/p900322/raw'\n",
- "fit_path = '/gpfs/exfel/data/scratch/mramilli/DATA/MID/p900322/gain_fit' # fit files directory\n",
+ "# fit_path = '/gpfs/exfel/data/scratch/mramilli/DATA/MID/p900322/gain_fit' # fit files directory\n",
"metadata_folder = '' # Directory containing calibration_metadata.yml when run by xfel-calibrate\n",
- "g0_run = 94 # Number of high gain\n",
+ "g0_runs = [94] # Number of high gain\n",
"g_map_old_dir = '/gpfs/exfel/data/user/mramilli/jungfrau/module_PSI_gainmaps/M302'\n",
"gain_map_file = '/gainMaps_M109_Burst_Fix_20230523.h5' # path\n",
"out_folder = \"/gpfs/exfel/data/scratch/ahmedk/test/jf_ff/gain_maps\"\n",
"karabo_id = 'MID_EXP_JF500K2' # karabo prefix of Jungfrau devices\n",
- "da_name = 'JNGFR02'\n",
+ "karabo_da = ['JNGFR02']\n",
"_fit_func = 'CHARGE_SHARING' # function used to fit the single photon peak\n",
"gains = [0, 1, 2]\n",
"# Parameter conditions\n",
@@ -40,8 +40,7 @@
"integration_time = -1 # the detector acquisition rate, use 0 to try to auto-determine\n",
"gain_setting = -1\n",
"gain_mode = -1\n",
- "memory_cells = 16\n",
- "mode = 'Burst'\n",
+ "memory_cells = -1\n",
"adc_fit = True\n",
"is_strixel = False\n",
"# Condition limits\n",
@@ -55,11 +54,11 @@
"db_output = False\n",
"send_bpix = False\n",
"gain_map_name = 'gain_map_g0'\n",
- "\n",
"g0_name_new = 'gainMap_fit' # name of the data structure in the fit files\n",
"E_ph = 8.048 # photon energy of the peak fitted\n",
"badpixel_threshold_sigma = 3. # number of std in gain distribution to mark bad pixels\n",
"_roi = [0, 1024, 0, 256] # ROI to consider to evaluate gain distribution (for bad pixels evaluation)\n",
+ "control_src_template = '{}/DET/CONTROL'\n",
"\n",
"# CALCAT API parameters\n",
"cal_db_interface = \"tcp://max-exfl-cal-001:8020\" # the database interface to use\n",
@@ -93,6 +92,8 @@
"from pathlib import Path\n",
"from XFELDetAna.plotting.histogram import histPlot\n",
"from XFELDetAna.plotting.heatmap import heatmapPlot\n",
+ "from cal_tools.jungfraulib import JungfrauCtrl\n",
+ "from extra_data import RunDirectory\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n"
]
@@ -110,14 +111,35 @@
"metadata": {},
"outputs": [],
"source": [
- "if memory_cells > 1:\n",
- " mode = 'Burst'\n",
- "\n",
- "fit_path = Path(fit_path)\n",
- "\n",
+ "g0_run = g0_runs[0]\n",
+ "fit_path = Path(out_folder) # TODO\n",
+ "da_name = karabo_da[0]\n",
"proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]\n",
"file_loc = f\"proposal:{proposal} runs:{g0_run}\"\n",
- "report = get_report(metadata_folder)"
+ "report = get_report(metadata_folder)\n",
+ "in_folder = Path(in_folder)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "run_folder = in_folder / f'r{g0_run:04d}'\n",
+ "control_src = control_src_template.format(karabo_id, karabo_da[0])\n",
+ "\n",
+ "ctrl_data = JungfrauCtrl(RunDirectory(run_folder), control_src)\n",
+ "\n",
+ "if memory_cells < 0:\n",
+ " memory_cells, _ = ctrl_data.get_memory_cells()\n",
+ "\n",
+ "print(f'Memory cells: {memory_cells}')\n",
+ "\n",
+ "if memory_cells > 1:\n",
+ " mode = 'Burst'\n",
+ "else:\n",
+ " mode = 'Single'"
]
},
{
@@ -135,7 +157,8 @@
" if is_strixel:\n",
" g0_new = np.transpose(from_strixel(np.transpose(g0_new)))\n",
"\n",
- " print(g0_new.shape)\n",
+ " print(\"Fit gain data shape:\", g0_new.shape)\n",
+ "\n",
" if integration_time < 0 and 'integration_time' in f.attrs.keys():\n",
" integration_time = np.float32(f.attrs['integration_time'])\n",
" print(f'Integration time: {integration_time}')\n",
@@ -147,7 +170,7 @@
" print(f'bias voltage: {bias_voltage}')\n",
" else:\n",
" print(f'bias voltage not found!\\nUsing default value {bias_voltage}')\n",
- " \n",
+ "\n",
" # if 'karabo_id' in f.attrs.keys():\n",
" # karabo_id = str(f.attrs['karabo_id'])\n",
" # print(f'Karabo ID: {karabo_id}')\n",
diff --git a/notebooks/Jungfrau/gainCal_JF_Create_Spectra_Histos.ipynb b/notebooks/Jungfrau/gainCal_JF_Create_Spectra_Histos.ipynb
index e9076a528c6c55d0f7b1ac91fb0025ed0bf1eccf..4dbe7b8aa4867430272e543529b48c4f0622fdf4 100644
--- a/notebooks/Jungfrau/gainCal_JF_Create_Spectra_Histos.ipynb
+++ b/notebooks/Jungfrau/gainCal_JF_Create_Spectra_Histos.ipynb
@@ -19,30 +19,56 @@
"metadata": {},
"outputs": [],
"source": [
- "instrument = 'HED'\n",
- "cycle = 202231\n",
- "proposal = 900297\n",
- "in_folder = f'/gpfs/exfel/exp/{instrument}/{cycle}/p{proposal:06d}/raw'\n",
- "fit_out_folder = f'/gpfs/exfel/data/scratch/ahmedk/test/FFDATA/{instrument}/p{proposal:06d}/gain_fit'\n",
+ "in_folder = '/gpfs/exfel/exp/HED/202231/p900297/raw'\n",
+ "out_folder = '/gpfs/exfel/data/scratch/ahmedk/test/FFDATA/HED/p900297/gain_fit'\n",
+ "metadata_folder = '' # Directory containing calibration_metadata.yml when run by xfel-calibrate\n",
"g0_runs = [26, ] # can be a list of runs\n",
"sensor_size = [512, 1024] # size of the array in the 'row' and 'col' dimensions\n",
"chunked_trains = 1000 # Number of trains per chunk.\n",
"gains = [0, 1, 2] # gain bit values\n",
- "mod_n = 1 # module number\n",
"karabo_id = 'HED_IA1_JF500K1' # karabo prefix of Jungfrau devices\n",
- "da_name = f'JNGFR{mod_n:02d}'\n",
+ "karabo_da = [\"JNGFR01\"]\n",
+ "da_template = 'JNGFR{:02d}'\n",
"creation_time = \"\" # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC e.g. \"2022-06-28 13:00:00\"\n",
- "bias_voltage = 180. # bias voltage - Set to -1 to derive from CONTROL file.\n",
- "integration_time = 10. # integration time - Set to -1 to derive from CONTROL file.\n",
- "gain_mode = 0 # number of memory cells - Set to -1 to derive from CONTROL file.\n",
- "memory_cells = 1 # number of memory cells - Set to -1 to derive from CONTROL file.\n",
- "sc_start = 15 # storage cell start value - should be derived from CONTROL file\n",
+ "bias_voltage = -1 # bias voltage - Set to -1 to derive from CONTROL file.\n",
+ "integration_time = -1 # integration time - Set to -1 to derive from CONTROL file.\n",
+ "gain_mode = -1 # number of memory cells - Set to -1 to derive from CONTROL file.\n",
+ "gain_setting = -1\n",
+ "\n",
+ "memory_cells = -1 # number of memory cells - Set to -1 to derive from CONTROL file.\n",
+ "sc_start = -1 # storage cell start value - should be derived from CONTROL file\n",
"h_bins_s = 0.2 # bin width in ADC units of the histogram\n",
"h_range = (-5., 10.) # histogram range in ADC units\n",
"block_size = [128, 64] # dimension of the chunks in 'row' and 'col'\n",
- "det_src = f'{karabo_id}/DET/{da_name}:daqOutput'\n",
- "control_src = f'{karabo_id}/DET/CONTROL'\n",
- "extra_dims = ['cell', 'row', 'col'] # labels for the DataArray dims after the first"
+ "det_src_template = '{}/DET/{}:daqOutput'\n",
+ "control_src_template = '{}/DET/CONTROL'\n",
+ "extra_dims = ['cell', 'row', 'col'] # labels for the DataArray dims after the first\n",
+ "fit_path = '/gpfs/exfel/data/scratch/mramilli/DATA/MID/p900322/gain_fit' # fit files directory\n",
+ "gain_map_file = '/gainMaps_M109_Burst_Fix_20230523.h5' # path\n",
+ "_fit_func = 'CHARGE_SHARING' # function used to fit the single photon peak\n",
+ "\n",
+ "mode = 'Burst' # TODO: CONFIRM MY DECISION TO REMOVE MODE DECLARATION FROM ALL NOTEBOOKS.\n",
+ "adc_fit = True\n",
+ "is_strixel = False\n",
+ "# Condition limits\n",
+ "bias_voltage_lims = [0, 200]\n",
+ "integration_time_lims = [0.1, 1000]\n",
+ "\n",
+ "g_map_old = '' # '/gpfs/exfel/data/user/mramilli/jungfrau/module_PSI_gainmaps/M302/gainMaps_M302_2022-02-23.h5' # old gain map file path to calculate gain ratios G0/G1 and G1/G2. Set to \"\" to get last gain constant from the database.\n",
+ "old_gain_dataset_name = 'gain_map_g0' # name of the data structure in the old gain map\n",
+ "correct_offset = True # correct the photon peak value with a pedestal fit position\n",
+ "db_output = False\n",
+ "send_bpix = False\n",
+ "gain_map_name = 'gain_map_g0'\n",
+ "\n",
+ "g0_name_new = 'gainMap_fit' # name of the data structure in the fit files\n",
+ "E_ph = 8.048 # photon energy of the peak fitted\n",
+ "badpixel_threshold_sigma = 3. # number of std in gain distribution to mark bad pixels\n",
+ "_roi = [0, 1024, 0, 256] # ROI to consider to evaluate gain distribution (for bad pixels evaluation)\n",
+ "\n",
+ "# CALCAT API parameters\n",
+ "cal_db_interface = \"tcp://max-exfl-cal-001:8020\" # the database interface to use\n",
+ "cal_db_timeout = 180000"
]
},
{
@@ -67,7 +93,20 @@
"from cal_tools.jungfrau import jungfrau_ff\n",
"from cal_tools.jungfraulib import JungfrauCtrl\n",
"from cal_tools.tools import calcat_creation_time\n",
- "from cal_tools.step_timing import StepTimer"
+ "from cal_tools.step_timing import StepTimer\n",
+ "from pathlib import Path"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "in_folder = Path(in_folder)\n",
+ "out_folder = Path(out_folder)\n",
+ "out_folder.mkdir(exist_ok=True, parents=True)\n",
+ "print(f'Creating directory {out_folder}')"
]
},
{
@@ -76,9 +115,9 @@
"metadata": {},
"outputs": [],
"source": [
- "if not os.path.exists(fit_out_folder):\n",
- " os.makedirs(fit_out_folder)\n",
- " print('Creating directory {}'.format(fit_out_folder))"
+ "det_src = det_src_template.format(karabo_id, karabo_da[0])\n",
+ "control_src = control_src_template.format(karabo_id, karabo_da[0])\n",
+ "da_name = karabo_da[0] # TODO fix this"
]
},
{
@@ -89,8 +128,8 @@
"source": [
"begin_stuff = time.localtime()\n",
"\n",
- "from pathlib import Path\n",
- "first_run_folder = Path(in_folder) / f'r{g0_runs[0]:04d}'\n",
+ "\n",
+ "first_run_folder = in_folder / f'r{g0_runs[0]:04d}'\n",
"ctrl_data = JungfrauCtrl(RunDirectory(first_run_folder), control_src)\n",
"\n",
"if memory_cells < 0 and sc_start < 0:\n",
@@ -121,7 +160,7 @@
"h_n_bins = int((h_range[1] - h_range[0])/h_bins_s)\n",
"h_bins = np.linspace(h_range[0], h_range[1], h_n_bins)\n",
"\n",
- "h_spectra = np.zeros((h_n_bins-1, memory_cells, sensor_size[0], sensor_size[1]), dtype=np.int)\n",
+ "h_spectra = np.zeros((h_n_bins-1, memory_cells, sensor_size[0], sensor_size[1]), dtype=np.int32)\n",
"edges = None"
]
},
@@ -131,93 +170,94 @@
"metadata": {},
"outputs": [],
"source": [
- "def read_conditions_and_retrieve_constants(\n",
- " run_dc,\n",
- " run_number,\n",
- " integration_time,\n",
- " bias_voltage,\n",
- " gain_mode,\n",
- " creation_time,\n",
- "):\n",
- " \"\"\"Retrieve offset and noise calibration constants\n",
- " for one jungfrau module.\n",
- " Args:\n",
- " run_dc(extra_data.DataCollection):\n",
- " run_number(int):\n",
- " integration_time(float)\n",
- " bias_voltage(float):\n",
- " gain_mode(int):\n",
- " creation_time():\n",
- " Returns:\n",
- " ndarray, ndarray: Offset and Noise calibration constants.\n",
- " \"\"\"\n",
- " \n",
- " ## Retrieve DB constants\n",
- " ctrl_data = JungfrauCtrl(run_dc, control_src)\n",
- " if integration_time < 0:\n",
- " integration_time = ctrl_data.get_integration_time()\n",
- " if bias_voltage < 0:\n",
- " bias_voltage = ctrl_data.get_bias_voltage()\n",
- " # gain_setting = ctrl_data.get_gain_setting() TODO: Where is the gain_setting?\n",
- " if gain_mode < 0:\n",
- " gain_mode = ctrl_data.get_gain_mode()\n",
- "\n",
- " # TODO: clarify this way in getting the control data.\n",
- " # We are using RUN. We can of course update and start using `as_single_value` method.\n",
- " print(f\"Memory cells: {memory_cells:d}\")\n",
- " print(f\"Exposure Time: {integration_time:1.2f} us\")\n",
- " print(f\"Bias Voltage: {bias_voltage:3.1f} V\")\n",
- "\n",
- " # Run's creation time:\n",
- " creation_time = calcat_creation_time(in_folder, run_number, creation_time)\n",
- " print(f\"Creation time: {creation_time}\")\n",
- "\n",
- " jf_cal = JUNGFRAU_CalibrationData(\n",
- " detector_name=karabo_id,\n",
- " sensor_bias_voltage=bias_voltage,\n",
- " event_at=creation_time,\n",
- " modules=da_name,\n",
- " memory_cells=memory_cells,\n",
- " integration_time=integration_time,\n",
- " # No gain setting 1? I use 0 for now as I don't see gain_setting used anywhere for the 3 notebooks.\n",
- " gain_setting=0,\n",
- " gain_mode=gain_mode,\n",
- " client=rest_cfg.calibration_client(),\n",
+ "# def read_conditions_and_retrieve_constants(\n",
+ "# run_dc,\n",
+ "# run_number,\n",
+ "# integration_time,\n",
+ "# bias_voltage,\n",
+ "# gain_mode,\n",
+ "# creation_time,\n",
+ "# ):\n",
+ "# \"\"\"Retrieve offset and noise calibration constants\n",
+ "# for one jungfrau module.\n",
+ "# Args:\n",
+ "# run_dc(extra_data.DataCollection):\n",
+ "# run_number(int):\n",
+ "# integration_time(float)\n",
+ "# bias_voltage(float):\n",
+ "# gain_mode(int):\n",
+ "# creation_time():\n",
+ "# Returns:\n",
+ "# ndarray, ndarray: Offset and Noise calibration constants.\n",
+ "# \"\"\"\n",
+ "run_folder = in_folder / f'r{g0_runs[0]:04d}' # first run for now\n",
+ "run_number = g0_runs[0]\n",
+ "## Retrieve DB constants\n",
+ "ctrl_data = JungfrauCtrl(RunDirectory(run_folder), control_src)\n",
+ "if integration_time < 0:\n",
+ " integration_time = ctrl_data.get_integration_time()\n",
+ "if bias_voltage < 0:\n",
+ " bias_voltage = ctrl_data.get_bias_voltage()\n",
+ "# gain_setting = ctrl_data.get_gain_setting() TODO: Where is the gain_setting?\n",
+ "if gain_mode < 0:\n",
+ " gain_mode = ctrl_data.get_gain_mode()\n",
+ "\n",
+ "# TODO: clarify this way in getting the control data.\n",
+ "# We are using RUN. We can of course update and start using `as_single_value` method.\n",
+ "print(f\"Memory cells: {memory_cells:d}\")\n",
+ "print(f\"Exposure Time: {integration_time:1.2f} us\")\n",
+ "print(f\"Bias Voltage: {bias_voltage:3.1f} V\")\n",
+ "\n",
+ "# Run's creation time:\n",
+ "creation_time = calcat_creation_time(in_folder, run_number, creation_time)\n",
+ "print(f\"Creation time: {creation_time}\")\n",
+ "\n",
+ "jf_cal = JUNGFRAU_CalibrationData(\n",
+ " detector_name=karabo_id,\n",
+ " sensor_bias_voltage=bias_voltage,\n",
+ " event_at=creation_time,\n",
+ " modules=karabo_da[0],\n",
+ " memory_cells=memory_cells,\n",
+ " integration_time=integration_time,\n",
+ " # No gain setting 1? I use 0 for now as I don't see gain_setting used anywhere for the 3 notebooks.\n",
+ " gain_setting=0,\n",
+ " gain_mode=gain_mode,\n",
+ " client=rest_cfg.calibration_client(),\n",
+ ")\n",
+ "jf_metadata = jf_cal.metadata(calibrations=[\"Offset10Hz\", \"Noise10Hz\"])\n",
+ "# TODO: display CCV timestamp\n",
+ "const_data = jf_cal.ndarray_map(metadata=jf_metadata)\n",
+ "\n",
+ "constants_data = const_data[karabo_da[0]]\n",
+ "\n",
+ "if \"Offset10Hz\" in constants_data:\n",
+ " offset_map = np.array(\n",
+ " np.transpose(np.swapaxes(constants_data[\"Offset10Hz\"], 0, 1)),\n",
+ " dtype=np.float32\n",
" )\n",
- " jf_metadata = jf_cal.metadata(calibrations=[\"Offset10Hz\", \"Noise10Hz\"])\n",
- " # TODO: display CCV timestamp\n",
- " const_data = jf_cal.ndarray_map(metadata=jf_metadata)\n",
- "\n",
- " constants_data = const_data[da_name]\n",
- "\n",
- " if \"Offset10Hz\" in constants_data:\n",
- " offset = np.array(\n",
- " np.transpose(np.swapaxes(constants_data[\"Offset10Hz\"], 0, 1)),\n",
- " dtype=np.float32\n",
- " )\n",
- " print(f'Retrieved Offset Map {offset.shape}')\n",
- " else:\n",
- " # Create offset empty constant.\n",
- " offset = np.zeros((3, memory_cells, 512, 1024), dtype=np.float32)\n",
- " print(\"Offset map not found\")\n",
- "\n",
- " if \"Noise10Hz\" in constants_data:\n",
- " noise = np.array(\n",
- " np.transpose(np.swapaxes(constants_data[\"Noise10Hz\"], 0, 1)),\n",
- " dtype=np.float32\n",
- " )\n",
- " print(f'Retrieved Noise Map {noise.shape}')\n",
- " else:\n",
- " noise = None\n",
- " print(\"Noise map not found\")\n",
- "\n",
- " if memory_cells > 1:\n",
- " # move from x, y, cell, gain to cell, x, y, gain\n",
- " offset = np.moveaxis(offset, [0, 1], [1, 2])\n",
- " else:\n",
- " offset = np.squeeze(offset)\n",
- "\n",
- " return offset, noise"
+ " print(f'Retrieved Offset Map {offset_map.shape}')\n",
+ "else:\n",
+ " # Create offset empty constant.\n",
+ " offset_map = np.zeros((3, memory_cells, 512, 1024), dtype=np.float32)\n",
+ " print(\"Offset map not found\")\n",
+ "\n",
+ "if \"Noise10Hz\" in constants_data:\n",
+ " noise_map = np.array(\n",
+ " np.transpose(np.swapaxes(constants_data[\"Noise10Hz\"], 0, 1)),\n",
+ " dtype=np.float32\n",
+ " )\n",
+ " print(f'Retrieved Noise Map {noise_map.shape}')\n",
+ "else:\n",
+ " noise_map = None\n",
+ " print(\"Noise map not found\")\n",
+ "\n",
+ "if memory_cells > 1:\n",
+ " # move from x, y, cell, gain to cell, x, y, gain\n",
+ " offset_map = np.moveaxis(offset_map, [0, 1], [1, 2])\n",
+ "else:\n",
+ " offset_map = np.squeeze(offset_map)\n",
+ "\n",
+ " # return offset, noise"
]
},
{
@@ -299,22 +339,22 @@
" \n",
" raw_fname_tmp = f'/RAW-R{g0_run:04d}-{da_name}' + '-S{:05d}.h5'\n",
"\n",
- " run_folder = in_folder + f'/r{g0_run:04d}'\n",
+ " run_folder = in_folder / f'r{g0_run:04d}'\n",
"\n",
- " run_dc = open_run(proposal=proposal, run=g0_run, data='raw')\n",
+ " run_dc = RunDirectory(run_folder)\n",
"\n",
" run_dc.info()\n",
"\n",
" run_dc = run_dc.select(selection)\n",
"\n",
- " offset_map, noise_map = read_conditions_and_retrieve_constants(\n",
- " run_dc=run_dc,\n",
- " run_number=g0_run,\n",
- " integration_time=integration_time,\n",
- " bias_voltage=bias_voltage,\n",
- " gain_mode=gain_mode,\n",
- " creation_time=creation_time,\n",
- " )\n",
+ " # offset_map, noise_map = read_conditions_and_retrieve_constants(\n",
+ " # run_dc=run_dc,\n",
+ " # run_number=g0_run,\n",
+ " # integration_time=integration_time,\n",
+ " # bias_voltage=bias_voltage,\n",
+ " # gain_mode=gain_mode,\n",
+ " # creation_time=creation_time,\n",
+ " # )\n",
"\n",
" ## Offset subtraction & Histogram filling\n",
" ### performs offset subtraction and fills the histogram\n",
@@ -390,7 +430,7 @@
"metadata": {},
"outputs": [],
"source": [
- "fout_h_path = f'{fit_out_folder}/{fout_temp}_Histo.h5'\n",
+ "fout_h_path = f'{out_folder}/{fout_temp}_Histo.h5'\n",
"hists = np.transpose(h_spectra) # this is to make it compatible with my other notebooks, can be removed\n",
"with h5file(fout_h_path, 'w') as fout_h: \n",
" print(f\"Saving histograms at {fout_h_path}.\")\n",
diff --git a/notebooks/Jungfrau/gainCal_JF_Fit_Spectra_Histos.ipynb b/notebooks/Jungfrau/gainCal_JF_Fit_Spectra_Histos.ipynb
index 367c0c6ac93a20d4b35b9333b7aa6a681a80ded6..a1299a58d9189e11576c35e973dcb7bbb61271e0 100644
--- a/notebooks/Jungfrau/gainCal_JF_Fit_Spectra_Histos.ipynb
+++ b/notebooks/Jungfrau/gainCal_JF_Fit_Spectra_Histos.ipynb
@@ -1,7 +1,6 @@
{
"cells": [
{
- "attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
@@ -18,33 +17,34 @@
"metadata": {},
"outputs": [],
"source": [
- "instrument = 'HED'\n",
- "cycle = 202231\n",
- "proposal = 900297\n",
- "g0_runs = [26,] # can be a list of runs\n",
- "dir_date_iso = '' ## string to save the date of the creation of the run \n",
- " ## to be used when injecting constants (third separate nb)\n",
- "filepath = f'/gpfs/exfel/exp/{instrument}/{cycle}/p{proposal:06d}/raw'\n",
+ "in_folder = \"/gpfs/exfel/exp/HED/202231/p900297/raw\"\n",
+ "out_folder = \"/gpfs/exfel/data/scratch/ahmedk/test/FFDATA/HED/p900297/gain_fit\"\n",
+ "runs = [26] # can be a list of runs\n",
+ "dir_date_iso = \"\" ## string to save the date of the creation of the run\n",
+ "\n",
+ "# to be used when injecting constants (third separate nb)\n",
"sensor_size = [512, 1024] # size of the array in the 'row' and 'col' dimensions\n",
"chunk_size = 10 # n of trains per chunk\n",
"block_size = [128, 64] # dimension of the chunks in 'row' and 'col'\n",
"mod_n = 1 # module number\n",
- "karabo_id = 'HED_IA1_JF500K1' # karabo prefix of Jungfrau devices\n",
- "da_name = f'JNGFR{mod_n:02d}'\n",
- "da_control = 'JNGFRCTRL00'\n",
+ "karabo_id = \"HED_IA1_JF500K1\" # karabo prefix of Jungfrau devices\n",
+ "da_name = f\"JNGFR{mod_n:02d}\"\n",
+ "da_control = \"JNGFRCTRL00\"\n",
+ "ctrl_source_template = \"{}/DET/CONTROL\" # template for control source name (filled with karabo_id_control)\n",
+ "karabo_id_control = \"\" # if control is on a different ID, set to empty string if it is the same a karabo-id\n",
"\n",
+ "bias_voltage = 180.0 # bias voltage - should be derived from CONTROL file\n",
+ "integration_time = 10.0 # integration time - should be derived from CONTROL file.\n",
"memory_cells = 1 # number of memory cells - should be derived from CONTROL file\n",
- "bias_voltage = 180. # bias voltage - should be derived from CONTROL file\n",
- "integration_time = 10. # integration time - should be derived from CONTROL file\n",
"sc_start = 15 # storage cell start value - should be derived from CONTROL file\n",
- "_fit_func = 'CHARGE_SHARING' ## which function will be used to fit the histogram\n",
- "_h_range = (200., 350.) # range of the histogram in x-axis units\n",
- "rebin = 1 \n",
+ "gain_mode = 0 # number of memory cells - Set to -1 to derive from CONTROL file.\n",
+ "\n",
+ "_fit_func = \"CHARGE_SHARING\" ## which function will be used to fit the histogram\n",
+ "_h_range = (200.0, 350.0) # range of the histogram in x-axis units\n",
+ "rebin = 1\n",
"# parameters for the peak finder\n",
- "n_sigma = 20. # n of sigma abov pedestal threshold\n",
- "ratio = 0.99 # ratio of the next peak amplitude in the peak_finder\n",
- "mode = 'Single'\n",
- "out_folder = f'/gpfs/exfel/data/scratch/mramilli/DATA/{instrument}/p{proposal:06d}/gain_fit'"
+ "n_sigma = 20.0 # n of sigma abov pedestal threshold\n",
+ "ratio = 0.99 # ratio of the next peak amplitude in the peak_finder"
]
},
{
@@ -53,16 +53,19 @@
"metadata": {},
"outputs": [],
"source": [
- "import numpy as np\n",
- "from h5py import File as h5file\n",
- "import logging\n",
+ "import multiprocessing\n",
"import time\n",
- "import os\n",
- "\n",
- "from pathlib import Path\n",
"from functools import partial\n",
- "import multiprocessing\n",
- "import jungfrau_ff"
+ "from logging import warning\n",
+ "from pathlib import Path\n",
+ "\n",
+ "from cal_tools.jungfrau import jungfrau_ff\n",
+ "import numpy as np\n",
+ "from cal_tools.jungfraulib import JungfrauCtrl\n",
+ "from extra_data import RunDirectory\n",
+ "from h5py import File as h5file\n",
+ "import pasha as psh\n",
+ "import time"
]
},
{
@@ -72,9 +75,11 @@
"outputs": [],
"source": [
"out_folder = Path(out_folder)\n",
- "if not os.path.exists(out_folder):\n",
- " os.makedirs(out_folder)\n",
- " print('Creating directory {}'.format(out_folder))"
+ "out_folder.mkdir(parents=True, exist_ok=True)\n",
+ "in_folder = Path(in_folder)\n",
+ "\n",
+ "if karabo_id_control == \"\":\n",
+ " karabo_id_control = karabo_id"
]
},
{
@@ -90,27 +95,55 @@
"metadata": {},
"outputs": [],
"source": [
+ "h_spectra = None\n",
+ "edges = None\n",
+ "noise_map = None\n",
+ "\n",
+ "ctrl_src = ctrl_source_template.format(karabo_id_control)\n",
+ "run_dc = RunDirectory(in_folder / f\"r{runs[0]:04d}\")\n",
+ "\n",
+ "ctrl_data = JungfrauCtrl(run_dc, ctrl_src)\n",
+ "\n",
+ "if memory_cells < 0:\n",
+ " memory_cells, sc_start = ctrl_data.get_memory_cells()\n",
+ "\n",
+ " mem_cells_name = \"single cell\" if memory_cells == 1 else \"burst\"\n",
+ " print(f\"Run is in {mem_cells_name} mode.\\nStorage cell start: {sc_start:02d}\")\n",
+ "else:\n",
+ " mem_cells_name = \"single cell\" if memory_cells == 1 else \"burst\"\n",
+ " print(\n",
+ " f\"Run is in manually set to {mem_cells_name} mode. With {memory_cells} memory cells\"\n",
+ " )\n",
+ "\n",
+ "mode = \"Single\"\n",
"if memory_cells > 1:\n",
- " mode = 'Burst'\n",
- "file_h_name = f'/R{g0_runs[0]:04d}_Gain_{mode}_Spectra_{da_name}_Histo.h5' # histogram file name\n",
+ " mode = \"Burst\"\n",
+ "\n",
+ "file_h_name = (\n",
+ " f\"R{runs[0]:04d}_Gain_{mode}_Spectra_{da_name}_Histo.h5\" # histogram file name\n",
+ ")\n",
"\n",
"begin_stuff = time.localtime()\n",
- "i_cut = file_h_name.find('_Histo')\n",
+ "i_cut = file_h_name.find(\"_Histo\")\n",
"fout_temp = file_h_name[:i_cut]\n",
+ "fout_temp += f\"_{_fit_func}_Fit\"\n",
+ "\n",
+ "print(f\"block_size: {block_size}\")\n",
+ "\n",
+ "if integration_time < 0:\n",
+ " integration_time = ctrl_data.get_integration_time()\n",
+ "if bias_voltage < 0:\n",
+ " bias_voltage = ctrl_data.get_bias_voltage()\n",
+ "# if gain_setting < 0:\n",
+ "# gain_setting = ctrl_data.get_gain_setting()\n",
+ "if gain_mode < 0:\n",
+ " gain_mode = ctrl_data.get_gain_mode()\n",
"\n",
- "fout_temp += f'_{_fit_func}_Fit'\n",
- "lout_n = fout_temp\n",
- "fout_n = fout_temp\n",
- "log_path = out_folder\n",
- "logger = logging.getLogger()\n",
- "fhandler = logging.FileHandler(filename='{}/{}.log'.format(log_path, lout_n), mode='a')\n",
- "formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')\n",
- "fhandler.setFormatter(formatter)\n",
- "logger.addHandler(fhandler)\n",
- "logger.setLevel(logging.DEBUG)\n",
- "\n",
- "logging.info('chunk_size: {}'.format(chunk_size))\n",
- "logging.info('block_size: {}'.format(block_size))"
+ "print(f\"Integration time is {integration_time} us\")\n",
+ "# print(f\"Gain setting is {gain_setting} (run settings: {ctrl_data.run_settings})\")\n",
+ "print(f\"Gain mode is {gain_mode} ({ctrl_data.run_mode})\")\n",
+ "print(f\"Bias voltage is {bias_voltage} V\")\n",
+ "print(f\"Number of memory cells are {memory_cells}\")"
]
},
{
@@ -126,50 +159,40 @@
"metadata": {},
"outputs": [],
"source": [
- "h_spectra = None\n",
- "edges = None\n",
- "noise_map = None\n",
- "# transposition here is just to make saved histo compatible with my other notebooks \n",
- "with h5file(out_folder + file_h_name, 'r') as f:\n",
- " logging.info(f'opening histo in file {file_h_name}')\n",
- " if 'histos' in f.keys():\n",
- " h_spectra = np.transpose(np.array(f['histos']))\n",
- " logging.info(f'histogram found: {h_spectra.shape}')\n",
+ "# transposition here is just to make saved histo compatible with my other notebooks\n",
+ "with h5file(out_folder / file_h_name, \"r\") as f:\n",
+ " print(f\"opening histo in file {file_h_name}\")\n",
+ " if \"histos\" in f.keys():\n",
+ " h_spectra = np.transpose(np.array(f[\"histos\"]))\n",
+ " print(f\"histogram found: {h_spectra.shape}\")\n",
" else:\n",
- " raise AttributeError('No histo in file!')\n",
- " \n",
- " edges = np.array(f['edges'])\n",
- " x = (edges[1:] + edges[:-1])/2.\n",
- " if 'noise_map' in f.keys():\n",
- " noise_map = np.array(f['noise_map'])\n",
- " logging.info(f'noise map found: {noise_map.shape}')\n",
+ " raise AttributeError(\"No histo in file!\")\n",
+ "\n",
+ " edges = np.array(f[\"edges\"])\n",
+ "\n",
+ " x = (edges[1:] + edges[:-1]) / 2.0\n",
+ " if \"noise_map\" in f.keys():\n",
+ " noise_map = np.array(f[\"noise_map\"])\n",
+ " print(f\"noise map found: {noise_map.shape}\")\n",
" else:\n",
- " logging.info('noise map not found!')\n",
- " \n",
- " if 'integration_time' in f.attrs.keys():\n",
- " integration_time = np.float32(f.attrs['integration_time'])\n",
+ " warning(\"noise map not found!\")\n",
+ "\n",
+ " print(\"Reading control data from histogram files.\")\n",
+ " if \"integration_time\" in f.attrs.keys():\n",
+ " integration_time = np.float32(f.attrs[\"integration_time\"])\n",
" else:\n",
- " logging.info(f'integration_time not found! using default value{integration_time}')\n",
- " if 'bias_voltage' in f.attrs.keys():\n",
- " bias_voltage = np.float32(f.attrs['bias_voltage'])\n",
+ " print(f\"integration_time is not found! using default value{integration_time}\")\n",
+ " if \"bias_voltage\" in f.attrs.keys():\n",
+ " bias_voltage = np.float32(f.attrs[\"bias_voltage\"])\n",
" else:\n",
- " logging.info(f'integration_time not found! using default value{integration_time}')\n",
- " \n",
- " if 'dir_date_iso' in f.attrs.keys():\n",
- " dir_date_iso = str(f.attrs['dir_date_iso'])\n",
+ " warning(f\"bias_voltage not found! using default value: {bias_voltage}\")\n",
+ "\n",
+ " if \"creation_time\" in f.attrs.keys():\n",
+ " dir_date_iso = str(f.attrs[\"creation_time\"])\n",
" if len(dir_date_iso) == 0:\n",
- " logging.info('warning: dir date is empty string')\n",
+ " warning(\"Dir date is empty string\")\n",
" else:\n",
- " logging.info('dir_date not found!')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "x, h_spectra = set_histo_range(x, h_spectra, _h_range) "
+ " print(\"dir_date not found!\")"
]
},
{
@@ -182,20 +205,27 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {
- "scrolled": false
- },
+ "metadata": {},
"outputs": [],
"source": [
- "chunks = jungfrau_ff.chunk_Multi([h_spectra, noise_map], block_size)\n",
- "print(np.array(chunks).shape)\n",
- "\n",
+ "chunks = jungfrau_ff.chunk_Multi([h_spectra], block_size)\n",
"pool = multiprocessing.Pool()\n",
"\n",
- "partial_fit = partial(jungfrau_ff.fit_histogram, x, _fit_func, n_sigma, rebin, ratio)\n",
+ "st = time.perf_counter()\n",
"\n",
- "logging.info('starting spectra fit')\n",
+ "partial_fit = partial(\n",
+ " jungfrau_ff.fit_histogram,\n",
+ " x,\n",
+ " _fit_func,\n",
+ " n_sigma,\n",
+ " rebin,\n",
+ " ratio,\n",
+ " noise_map,\n",
+ ")\n",
+ "\n",
+ "print(\"starting spectra fit\")\n",
"r_maps = pool.map(partial_fit, chunks)\n",
+ "print(\"r_maps calculation\", time.perf_counter() - st)\n",
"\n",
"g0_map = np.zeros((memory_cells, sensor_size[0], sensor_size[1]), dtype=np.float32)\n",
"sigma_map = np.zeros((memory_cells, sensor_size[0], sensor_size[1]), dtype=np.float32)\n",
@@ -204,26 +234,31 @@
"\n",
"for i, r in enumerate(r_maps):\n",
" g0_chk, sigma_chk, chi2ndf_chk, alpha_chk = r\n",
- " \n",
- " n_blocks_col = int(g0_map.shape[-1]/block_size[1])\n",
- " irow = int(np.floor(i/n_blocks_col)) * block_size[0]\n",
- " icol = i%n_blocks_col * block_size[1]\n",
- " \n",
- " g0_map[..., irow:irow+block_size[0], icol:icol+block_size[1]] = g0_chk\n",
- " sigma_map[..., irow:irow+block_size[0], icol:icol+block_size[1]] = sigma_chk\n",
- " chi2ndf_map[..., irow:irow+block_size[0], icol:icol+block_size[1]] = chi2ndf_chk\n",
- " alpha_map[..., irow:irow+block_size[0], icol:icol+block_size[1]] = alpha_chk\n",
- " \n",
+ "\n",
+ " n_blocks_col = int(g0_map.shape[-1] / block_size[1])\n",
+ " irow = int(np.floor(i / n_blocks_col)) * block_size[0]\n",
+ " icol = i % n_blocks_col * block_size[1]\n",
+ "\n",
+ " g0_map[..., irow : irow + block_size[0], icol : icol + block_size[1]] = g0_chk\n",
+ " sigma_map[..., irow : irow + block_size[0], icol : icol + block_size[1]] = sigma_chk\n",
+ " chi2ndf_map[\n",
+ " ..., irow : irow + block_size[0], icol : icol + block_size[1]\n",
+ " ] = chi2ndf_chk\n",
+ " alpha_map[..., irow : irow + block_size[0], icol : icol + block_size[1]] = alpha_chk\n",
+ "\n",
+ "print(\"loading r_maps calculation results\", time.perf_counter() - st)\n",
+ "\n",
"pool.close()\n",
- "logging.info('... done')"
+ "\n",
+ "print(\"... done\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## final steps ##\n",
- "### saving fit results ###"
+ "## Final steps\n",
+ "### Saving fit results"
]
},
{
@@ -232,33 +267,33 @@
"metadata": {},
"outputs": [],
"source": [
+ "fout_path = f\"{out_folder}/{fout_temp}.h5\"\n",
"\n",
- "fout_path = f'{out_folder}/{fout_n}.h5'\n",
- "logging.info(fout_path)\n",
- "with h5file(fout_path, 'w') as fout:\n",
- " logging.info('saving noise map ...')\n",
+ "with h5file(fout_path, \"w\") as fout:\n",
+ " print(\"saving noise map ...\")\n",
"\n",
" ##trasposition is to make it compatible with my other nb\n",
" dset_noi = fout.create_dataset(\"noise_map\", data=np.transpose(noise_map))\n",
- " logging.info('saving fit results ...')\n",
+ " print(\"saving fit results ...\")\n",
"\n",
" dset_chi2 = fout.create_dataset(\"chi2map\", data=np.transpose(chi2ndf_map))\n",
" dset_gmap_fit = fout.create_dataset(\"gainMap_fit\", data=np.transpose(g0_map))\n",
" dset_std = fout.create_dataset(\"sigmamap\", data=np.transpose(sigma_map))\n",
" dset_alpha = fout.create_dataset(\"alphamap\", data=np.transpose(alpha_map))\n",
+ " fout.attrs[\"memory_cells\"] = memory_cells # TODO: Why memory cells are not saved here. What about the other conditions??\n",
" fout.attrs[\"integration_time\"] = integration_time\n",
" fout.attrs[\"bias_voltage\"] = bias_voltage\n",
" fout.attrs[\"dir_date_iso\"] = dir_date_iso\n",
- " fout.attrs['karabo_id'] = karabo_id\n",
- " fout.attrs['da_name'] = da_name\n",
+ " fout.attrs[\"karabo_id\"] = karabo_id\n",
+ " fout.attrs[\"da_name\"] = da_name\n",
"\n",
- "logging.info('closing')"
+ "print(\"closing\")"
]
}
],
"metadata": {
"kernelspec": {
- "display_name": "Python 3",
+ "display_name": ".cal2_venv",
"language": "python",
"name": "python3"
},
@@ -276,5 +311,5 @@
}
},
"nbformat": 4,
- "nbformat_minor": 1
+ "nbformat_minor": 4
}
diff --git a/notebooks/Jungfrau/gainCal_JF_Fit_sendDB_NBC.ipynb b/notebooks/Jungfrau/gainCal_JF_Fit_sendDB_NBC.ipynb
index 9d7d7d4975ac91e2fd7a5edafdd65d96875f33f9..f145000e831d458ec8f005e6d92ab94dbef4d95d 100755
--- a/notebooks/Jungfrau/gainCal_JF_Fit_sendDB_NBC.ipynb
+++ b/notebooks/Jungfrau/gainCal_JF_Fit_sendDB_NBC.ipynb
@@ -19,11 +19,11 @@
"out_folder = \"/gpfs/exfel/data/scratch/ahmedk/test/jf_ff/gain_maps\"\n",
"metadata_folder = '' # Directory containing calibration_metadata.yml when run by xfel-calibrate\n",
"\n",
- "g0_run = 94\n",
+ "g0_runs = [94]\n",
"\n",
"# Detector module parameters.\n",
"karabo_id = 'MID_EXP_JF500K2'\n",
- "da_name = 'JNGFR02'\n",
+ "karabo_da = ['JNGFR02']\n",
"\n",
"# Parameter conditions\n",
"bias_voltage = -1 # detector bias voltage\n",
@@ -85,12 +85,12 @@
"source": [
"if gain_setting == 1:\n",
" gain_map_name = 'gain_map_hg0'\n",
- "\n",
+ "da_name = karabo_da[0]\n",
"g_map_old_dir = Path(g_map_old_dir)\n",
"in_folder = Path(in_folder)\n",
"out_folder = Path(out_folder)\n",
"gain_file = g_map_old_dir / gain_map_file\n",
- "\n",
+ "g0_run = g0_runs[0]\n",
"\n",
"# Run's creation time:\n",
"creation_time = calcat_creation_time(in_folder, g0_run, creation_time)\n",
diff --git a/notebooks/Jungfrau/jungfrau_ff.py b/notebooks/Jungfrau/jungfrau_ff.py
deleted file mode 100644
index b152d9521291768fedee8ba9af5e2c1aafaa9649..0000000000000000000000000000000000000000
--- a/notebooks/Jungfrau/jungfrau_ff.py
+++ /dev/null
@@ -1,508 +0,0 @@
-import numpy as np
-from XFELDetAna.detectors.jungfrau.fitFuncs import (
- charge_sharing,
- chi_square_fit,
- double_peak,
- gauss,
- histogram_errors,
-)
-
-
-def chunk_trains(data, train_chk):
- """
- chunks data along the train dimension
-
- Args:
- * data (list): input data arrays, each at least 1-D
- * train_chk (int): number of trains per chunk
-
- Returns: list of chunked data
- """
- n_trains = int(data[0].shape[0])
- chunks = []
- for i in range(0, n_trains, train_chk):
- this_chunk = []
- this_chunk = [data[0][i:min(i+train_chk, n_trains)]]
- for j in range(1, len(data)):
- if data[j] is not None:
- this_chunk.append(data[j][i:min(i+train_chk, n_trains)])
- else:
- this_chunk.append(None)
-
- chunks.append(this_chunk)
- return chunks
-
-
-def chunk_Multi(data, block_size):
- """
- chunks input array along the 'row' and 'col' dimensions
-
- Args:
- * data (list): input data arrays, each at least 2-D
- * block_size (list, int): [chunk_row, chunk_col],
- dimension of the 2-D chunk
-
- Returns: list of chunked data
- """
- intervals_r = range(0, data[0].shape[-2], block_size[0])
- intervals_c = range(0, data[0].shape[-1], block_size[1])
-
- chunks = []
-
- for irow in intervals_r:
- for icol in intervals_c:
- this_chunk = []
- this_chunk = [
- data[0][
- ..., irow:irow + block_size[0],
- icol:icol + block_size[1]
- ]
- ]
- for j in range(1, len(data)):
- if data[j] is not None:
- this_chunk.append(
- data[j][..., irow:irow + block_size[0],
- icol:icol + block_size[1]]
- )
- else:
- this_chunk.append(None)
- chunks.append(this_chunk)
-
- return chunks
-
-
-def subtract_offset(offset_map, _inp):
- """
- perform offset subtraction on raw data
-
- Args:
- * offset_map (xarray, float): map with offset constants,
- with shape (3, memory_cells, row, col)
- * _inp (list): input data as:
- * _inp[0]: raw images, with shape (trains, memory_cells, row, col)
- * _inp[1]: gain bit map, with shape (trains, memory_cells, row, col)
-
- Returns: offset subtracted images
- """
- imgs = _inp[0]
- gbit = _inp[1]
-
- n_cells = imgs.shape[1]
-
- for cell in range(n_cells):
- this_cell_gbit = gbit.sel(cell=cell)
-
- this_cell_off = offset_map[:, cell]
-
- _o = np.choose(
- this_cell_gbit, (
- np.expand_dims(this_cell_off[0], axis=0),
- np.expand_dims(this_cell_off[1], axis=0),
- np.expand_dims(this_cell_off[2], axis=0)
- )
- )
-
- imgs.loc[dict(cell=cell)] -= _o
-
- return imgs
-
-
-def fill_histogram(h_bins, _inp):
- """
- Fills an histogram with shape (n_bins-1, memory_cells, n_rows, n_cols)
- from input images.
-
- Args:
- * h_bins (list, float): the binning of the x-axis
- * _inp (list): data to be histogrammed, as:
- * _inp[0]: images, with shape (trains, memory_cells, row, col)
-
- Returns: histogram bin counts, bin edges.
-
- """
- imgs = _inp[0]
- n_cells = imgs.shape[1]
- n_rows = imgs.shape[2]
- n_cols = imgs.shape[3]
-
- h_chk = np.zeros((len(h_bins)-1, n_cells, n_rows, n_cols), dtype=np.int)
-
- e_chk = None
-
- for cell in range(n_cells):
- for row in range(n_rows):
- for col in range(n_cols):
-
- this_cell = imgs.sel(cell=cell)
- this_pix = this_cell.isel(row=row, col=col)
-
- _h, _e = np.histogram(this_pix, bins=h_bins)
- h_chk[..., cell, row, col] += _h
-
- if e_chk is None:
- e_chk = np.array(_e)
-
- return h_chk, e_chk
-
-
-# peak finder to find the first photon peak and/or the pedestal peak
-# can/should be replaced with smthg more efficient
-def peak_position(x, h, thr=75, dist=30., ratio=0.4):
- """
- finds peaks in the histogram
- very rough, self-made function, should be replaced by better one
-
- Args:
- * x (list, float): histogram bin centers
- * h (list, float): histogram bin counts
- * thr (int): lower threshold along x-axis in ADC units to look for a peak
- * dist (float): minimal distance between the last found peak and the next
- * ratio (float): minimal ratio between the last found peak and the next
-
- Returns: list of sorted peak positions along the x-axis,
- index of the peak with lowes value along x-axis
- """
- imin = 0
- peak_bins = []
- if len(x) == len(h):
- for i in range(0, len(x)):
- if x[i] > thr:
- imin = i
- break
- sorted_i = np.argsort(h[imin:])
- peak_bins.append(sorted_i[-1] + imin)
- low_h = ratio * h[sorted_i[-1] + imin]
- for j in range(1, len(sorted_i)):
- ibin = sorted_i[-1 - j] + imin
- if h[ibin] > low_h:
- min_dist = (x[peak_bins] - x[ibin])**2. > dist**2.
- if np.all(min_dist):
- peak_bins.append(ibin)
- else:
- raise AttributeError(
- "axis and histogram arrays have "
- f"different length: {len(x)} and {len(h)}")
- return peak_bins, imin
-
-
-def fit_charge_sharing(x, y, yerr, sigma0, n_sigma, ratio):
- """
- fits histogram with single photon function with charge sharing tail
-
- Args:
- * x (array, float): x values
- * y (array, float): y values
- * yerr (array, float): errors of the y values
- * sigma0 (float): rough estimate of peak variance
- * n_sigma (int): to calculate threshold of the peak finder as
- thr = n_sigma * sigma0
- * ratio (float): ratio parameter of the peak finder
-
- Returns: peak position, peak variance, chi2/ndf,
- charge sharing parameter.
- """
-
- norm = np.sum(y) * (x[1] - x[0])
- alpha0 = 0.3
- thr = n_sigma * sigma0
- _peaks, _ = peak_position(x, y, thr, ratio=ratio)
-
- q0 = -1.
- if len(_peaks) > 0:
- q0 = np.ma.min(x[_peaks])
-
- q = -1.
- sigma = -1.
- alpha = -1.
- chi2ndf = -1.
-
- if q0 > -1.:
- limit_q = (q0-sigma0, q0+sigma0)
- limit_sigma = (0.1*sigma0, 2.*sigma0)
- limit_alpha = (0.01*alpha0, 1.)
- limit_norm = (np.sqrt(norm), 3.*norm)
- res = chi_square_fit(
- charge_sharing, x, y, yerr,
- fit_range=(0.5*q0, q0+1.5*sigma0),
- ped=0.,
- q=q0,
- sigma=sigma0,
- alpha=alpha0,
- norm=norm,
- limit_q=limit_q,
- limit_sigma=limit_sigma,
- limit_alpha=limit_alpha,
- limit_norm=limit_norm,
- fix_ped=True,
- print_level=0,
- pedantic=False
- )
- values = res[0]
- chi2 = res[2]
- ndf = res[3]
- is_valid = res[4]
-
- if is_valid:
- q = values['q']
- sigma = values['sigma']
- alpha = values['alpha']
- chi2ndf = chi2/ndf
- else:
- q = -1000.
- sigma = -1000.
- alpha = -1000.
- chi2ndf = -1000.
- return q, sigma, chi2ndf, alpha
-
-
-def fit_double_charge_sharing(x, y, yerr, sigma0, n_sigma, ratio):
- """
- fits histogram with sum of two single photon functions
- with charge sharing tail.
-
- Args:
- * x (array, float): x values
- * y (array, float): y values
- * yerr (array, float): errors of the y values
- * sigma0 (float): rough estimate of peak variance
- * n_sigma (int): to calculate threshold of the peak finder as
- thr = n_sigma * sigma0
- * ratio (float): ratio parameter of the peak finder
-
- Returns: peak position, peak variance, chi2/ndf,
- charge sharing parameter (all of them of the first peak)
-
- """
- norm = np.sum(y) * (x[1] - x[0])
- alpha0 = 0.3
- thr = n_sigma * sigma0
- _peaks, _ = peak_position(x, y, thr, ratio=ratio)
-
- qa = -1.
- if len(_peaks) > 0:
- qa = np.ma.min(x[_peaks])
-
- q = -1.
- sigma = -1.
- alpha = -1.
- chi2ndf = -1.
-
- if qa > -1.:
-
- qb = 1.1 * qa
- limit_q_a = (qa - sigma0, qa + sigma0)
- limit_sigma_a = (0.1 * sigma0, 2. * sigma0)
- limit_alpha = (0.01, 1.)
- limit_norm_a = (np.ma.sqrt(0.1 * norm), 2. * norm)
- limit_q_b = (qb - sigma0, qb + sigma0)
- limit_sigma_b = (0.1 * sigma0, 2. * sigma0)
- limit_norm_b = (np.ma.sqrt(0.1 * norm), 2 * norm)
-
- res = chi_square_fit(
- double_peak, x, y, yerr,
- fit_range=(0.5*qa, qb + sigma0),
- ped=0.,
- q_a=qa,
- sigma_a=sigma0,
- alpha=alpha0,
- norm_a=norm,
- q_b=qb,
- sigma_b=sigma0,
- norm_b=norm,
- limit_q_a=limit_q_a,
- limit_sigma_a=limit_sigma_a,
- limit_alpha=limit_alpha,
- limit_norm_a=limit_norm_a,
- limit_q_b=limit_q_b,
- limit_sigma_b=limit_sigma_b,
- limit_norm_b=limit_norm_b,
- fix_ped=True,
- print_level=0,
- pedantic=False
- )
- values = res[0]
- chi2 = res[2]
- ndf = res[3]
- is_valid = res[4]
-
- if is_valid:
- q = values['q_a']
- sigma = values['sigma_a']
- alpha = values['alpha']
- chi2ndf = chi2/ndf
- else:
- q = -1000.
- sigma = -1000.
- alpha = -1000.
- chi2ndf = -1000.
-
- return q, sigma, chi2ndf, alpha
-
-
-def fit_gauss(x, y, yerr, sigma0, n_sigma, ratio):
- """
- fits histogram with a gaussian function
-
- Args:
- * x (array, float): x values
- * y (array, float): y values
- * yerr (array, float): errors of the y values
- * sigma0 (float): rough estimate of peak variance
- * n_sigma (int): to calculate threshold of the peak finder
- as thr = n_sigma * sigma0
- * ratio (float): ratio parameter of the peak finder
-
- Returns: peak position, peak variance, chi2/ndf,
- charge sharing parameter (last one alway == 0)
- all of them are kept for compatibility with the wrap around function.
-
- """
- norm = np.sum(y) * (x[1] - x[0])/np.sqrt(2.*np.pi*sigma0**2)
- thr = n_sigma * sigma0
- _peaks, _ = peak_position(x, y, thr, ratio=ratio)
-
- q0 = -1000.
- if len(_peaks) > 0:
- q0 = np.ma.min(x[_peaks])
-
- q = -1.
- sigma = -1.
- alpha = -1.
- chi2ndf = -1.
-
- if q0 > -1000.:
- limit_amp = (0.1*norm, 10.*norm)
- limit_mean = (q0-sigma0, q0+sigma0)
- limit_sigma = (0.1*sigma0, 10.*sigma0)
-
- res = chi_square_fit(gauss, x, y, yerr,
- fit_range=[q0 - 2.*sigma0, q0 + 2.*sigma0],
- amp=norm,
- mean=q0,
- sigma=sigma0,
- limit_amp=limit_amp,
- limit_mean=limit_mean,
- limit_sigma=limit_sigma,
- print_level=0,
- pedantic=False)
- values = res[0]
- chi2 = res[2]
- ndf = res[3]
- is_valid = res[4]
-
- if is_valid:
- q = values['mean']
- sigma = values['sigma']
- alpha = 0.
- chi2ndf = chi2/ndf
- else:
- q = -1000.
- sigma = -1000.
- alpha = -1000.
- chi2ndf = -1000.
-
- return q, sigma, chi2ndf, alpha
-
-
-def rebin_histo(h, x, r):
-
- if len(x) % r == 0:
- h_out = None
- x_out = x[::r]
-
- for i in range(0, len(h), r):
- if h_out is not None:
- h_out = np.append(h_out, np.sum(h[i:i+r], axis=0))
- else:
- h_out = np.array(np.sum(h[i:i+r], axis=0))
-
- else:
- raise AttributeError('Rebin factor is not exact divisor of bins!')
-
- return h_out, x_out
-
-
-def set_histo_range(_x, _h, _h_range):
- """
- reduces the histogram bin axis to the range specified
-
- Args:
- * _x (array, float): the bin centers array
- * _h (array, integers): the histogram with shape
- (bins, cells, columns, row)
- * _h_range (tuple, float): the (liminf, limsup) of the desired range
-
- Returns: the new bin centers array and the new histogram
-
- """
-
- i_min = (np.abs(_x - _h_range[0])).argmin()
- i_max = (np.abs(_x - _h_range[1])).argmin()
-
- return _x[i_min:i_max], _h[i_min:i_max]
-
-
-def fit_histogram(x, _fit_func, n_sigma, rebin, ratio, _inp):
- """
- wrap around function for fitting of histogram
-
- Args:
- * x (list, float): - bin centers along x
- * _fit_func (string): which function to use for fit
- * CHARGE_SHARING: single peak with charge sharing tail
- * CHARGE_SHARING_2: sum of two CHARGE_SHARING
- * GAUSS: gaussian function
- * n_sigma (int): to calculate threshold of the peak finder
- as thr = n_sigma * sigma0
- * ratio (float): ratio parameter of the peak finder
- * _input (list): contains the data to preform the fit
- * _input[0]: histogram bin counts with shape
- (n_bins, memory_cells, row, col)
- * _input[1]: noise map with shape (3, memory_cells, row, col)
-
- Returns: map of peak values, map of peak variances,
- map of chi2/ndf, map of charge sharing parameter values
- """
-
- _funcs = dict(
- CHARGE_SHARING=fit_charge_sharing,
- GAUSS=fit_gauss,
- CHARGE_SHARING_2=fit_double_charge_sharing,
- )
-
- fit_func = _funcs[_fit_func]
-
- histo = _inp[0]
- noise = _inp[1]
-
- n_cells = histo.shape[1]
- n_rows = histo.shape[2]
- n_cols = histo.shape[3]
-
- sigma0 = 15.
-
- g0_chk = np.zeros((n_cells, n_rows, n_cols), dtype=np.float32)
- sigma_chk = np.zeros((n_cells, n_rows, n_cols), dtype=np.float32)
- chi2ndf_chk = np.zeros((n_cells, n_rows, n_cols), dtype=np.float32)
- alpha_chk = np.zeros((n_cells, n_rows, n_cols), dtype=np.float32)
-
- for cell in range(n_cells):
- for row in range(n_rows):
- for col in range(n_cols):
-
- y = histo[..., cell, row, col]
- y, _x = rebin_histo(y, x, rebin)
- yerr = histogram_errors(y)
-
- if noise[0, cell, row, col] > 0.:
- sigma0 = noise[0, cell, row, col]
-
- q, sigma, chi2ndf, alpha = fit_func(
- _x, y, yerr, sigma0, n_sigma, ratio)
-
- g0_chk[cell, row, col] = q
- sigma_chk[cell, row, col] = sigma
- chi2ndf_chk[cell, row, col] = chi2ndf
- alpha_chk[cell, row, col] = alpha
-
- return g0_chk, sigma_chk, chi2ndf_chk, alpha_chk
diff --git a/src/xfel_calibrate/calibrate.py b/src/xfel_calibrate/calibrate.py
index fa14605bf1068aa43b0e512332b607ed0afb0b74..6d2d3104921b1dc3e881b2807d3104bd58eeace6 100755
--- a/src/xfel_calibrate/calibrate.py
+++ b/src/xfel_calibrate/calibrate.py
@@ -576,12 +576,9 @@ def run(argv=None):
# Ensure files are opened as UTF-8 by default, regardless of environment.
if "readthedocs.org" not in sys.executable:
locale.setlocale(locale.LC_CTYPE, ('en_US', 'UTF-8'))
-
if argv is None:
argv = sys.argv
-
args, nb_details = parse_argv_and_load_nb(argv)
-
concurrency = nb_details.concurrency
concurrency_par = args["concurrency_par"] or concurrency['parameter']
if concurrency_par == concurrency['parameter']: