diff --git a/cal_tools/cal_tools/agipdutils_ff.py b/cal_tools/cal_tools/agipdutils_ff.py
deleted file mode 100644
index 379b86f270c8bfbb9dc265b44cc471c4dc2e58ab..0000000000000000000000000000000000000000
--- a/cal_tools/cal_tools/agipdutils_ff.py
+++ /dev/null
@@ -1,243 +0,0 @@
-from iminuit import Minuit
-import numpy as np
-
-from cal_tools.enums import BadPixelsFF
-
-
-def all_in(msk, bits):
- return msk.astype(np.uint) & bits == bits
-
-
-def any_in(msk, bits):
- return msk.astype(np.uint) & bits > 0
-
-
-def gauss_fun(x, norm=1, mean=0, sigma=1):
- """
- Return value of Gaussian function
-
- :param x: Argument (float of 1D array) of Gaussian function
- :param norm: Normalization of Gaussian function
- :param mean: Mean parameter
- :param sigma: Sigma parameter
- :return: Value of gaussian function.
- """
- return norm * np.exp(-1 / 2 * ((x - mean) / sigma) ** 2) / (
- sigma * np.sqrt(2 * np.pi))
-
-
-def multi_gauss(x, ng=4, *p):
- """
- Sum of ng Gaussian functions
- :param x: Argument (float of 1D array) of the function
- :param ng: Number of Gaussian functions
- :param p: List of parameters (norm1,mean1,sigma1,norm2,mean2,sigma2,...)
- """
- r = 0.
- for i in range(ng):
- r += gauss_fun(x, *p[i * 3:(i + 1) * 3])
- return r
-
-
-def get_n_m_s(x, y, x_range):
- """
- Return statistical parameters of selected part of a histogram.
- ToDo: Check if wq.median works better than mean
-
- :param x: Center of bins of the histogram
- :param y: Value of bins of the histogram
- :param x_range: x range to be considered
- :return: Sum of histogram, Mean value, Standard deviation,
- List of selected bins
- """
- sel = (x >= x_range[0]) & (x < x_range[1])
- h_sum = np.sum(y[sel])
- h_norm = y[sel] / h_sum
- h_mean = np.sum(h_norm * x[sel])
- # h_med = wq.median(x[sel], y[sel])
- h_sqr = (x[sel] - h_mean) ** 2
- h_std = np.sqrt(np.sum(h_norm * h_sqr))
-
- return h_sum, h_mean, h_std, sel
-
-
-def get_starting_par(xe, ye, limits,
- n_peaks=3, f_lim=2):
- """
- Estimate starting parameters for Gaussian fit of several peaks.
-
- :param xe: Center of bins of the histogram
- :param ye: Value of bins of the histogram
- :param limits: Position of each peak ((left1, right1),
- (left2, right2), ...) to be considered.
- :param n_peaks: Number of peaks
- :param f_lim: Limits in units of standard deviation to consider
- """
- pars = {}
- shapes = []
- for peak in range(n_peaks):
- n, m, rms, idx = get_n_m_s(xe, ye, limits[peak])
- limits2 = [m - f_lim * rms, m + f_lim * rms]
- n, m, rms, idx = get_n_m_s(xe, ye, limits2)
- shapes.append((n, m, rms, idx))
-
- pars.update({f'g{peak}sigma': rms,
- f'g{peak}n': float(n),
- f'g{peak}mean': m
- }
- )
- return pars, shapes
-
-
-def fit_n_peaks(x, y, pars, x_range,
- do_minos=False,
- n_peaks=4,
- fix_d01=True
- ):
- """
- Fit histogram with n-Gaussian function.
-
- :param x: Center of bins of the histogram
- :param y: Value of bins of the histogram
- :param pars: Dictionary of initial parameters for fitting
- :param x_range: x Range to be considered for the fitting
- :param do_minos: Run Minos if True
- :param n_peaks: Number of Gaussian peaks to fit (min 2, max 4)
- :param fix_d01: Fix position of peaks to the distance between noise and
- first photon peak.
- :return: minuit object
- """
- sel = (x >= x_range[0]) & (x < x_range[1])
-
- # Square of bin errors
- yrr2 = np.copy(y[sel])
- yrr2[yrr2 == 0] = 1 # bins with zero events have error=1
-
- if fix_d01:
- pars['fix_g2mean'] = True
- pars['fix_g3mean'] = True
-
- if n_peaks < 4:
- pars['g3n'] = 0
- pars['fix_g3n'] = True
- pars['g3sigma'] = 1
- pars['fix_g3sigma'] = True
- pars['fix_g3mean'] = True
-
- if n_peaks < 3:
- pars['g2n'] = 0
- pars['fix_g2n'] = True
- pars['g2sigma'] = 1
- pars['fix_g2sigma'] = True
- pars['fix_g2mean'] = True
-
- def chi2_f(g0n, g0mean, g0sigma,
- g1n, g1mean, g1sigma,
- g2n, g2mean, g2sigma,
- g3n, g3mean, g3sigma, ):
-
- d01 = (g1mean - g0mean)
-
- if 'fix_g2mean' in pars and pars['fix_g2mean']:
- g2mean = g0mean + d01 * 2
-
- if 'fix_g3mean' in pars and pars['fix_g3mean']:
- g3mean = g0mean + d01 * 3
-
- if g3n == 0:
- n_peaks = 3
- elif g2n == 0:
- n_peaks = 2
- else:
- n_peaks = 4
-
- yt = multi_gauss(x[sel], n_peaks,
- g0n, g0mean, g0sigma,
- g1n, g1mean, g1sigma,
- g2n, g2mean, g2sigma,
- g3n, g3mean, g3sigma)
- return np.nansum((yt - y[sel]) ** 2 / yrr2)
-
- minuit = Minuit(chi2_f, **pars, pedantic=False)
- minuit.migrad()
-
- if do_minos:
- if minuit.get_fmin().is_valid:
- minuit.minos()
-
- return minuit
-
-
-def set_par_limits(pars, peak_range, peak_norm_range,
- peak_width_range, n_peaks=4):
- """
- Set limits on initial fit parameters based on given values
-
- ToDo: Check if limits on peak position should be given for all peaks
-
- :param pars: Dictionary of initial fit parameters
- :param peak_range: Limits on peak positions
- :param peak_norm_range: Limits on normalization of Gaussian peaks
- :param peak_width_range: Limits on width of Gaussian peaks
- :param n_peaks: Number of Gaussian peaks
- """
- pars.update({f'limit_g0mean': peak_range[0]})
- for peak in range(n_peaks):
- pars.update({f'limit_g{peak}n': peak_norm_range[peak],
- f'limit_g{peak}sigma': peak_width_range[peak],
- })
-
-
-def is_inside(x, a, b):
- """
- Check if float x is inside of range (a,b)
- """
- return (x >= a) & (x <= b)
-
-
-def get_mask(fit_summary, peak_lim, d0_lim, chi2_lim, peak_width_lim):
- """
- Calculate Bad pixels mask based on fit results and given limits
-
- :param fit_summary: Dictionary of the fit output from Multi-Gaussian fit
- :param peak_lim: Limits on noise peak position
- :param d0_lim: Limits on distance between noise and first photon peak
- :param chi2_lim: Limits on reduced chi^2 value
- :param peak_width_lim: Limits on noise peak width
- :return: Bad pixel mask
- """
- if not fit_summary['is_valid']:
- return BadPixelsFF.FIT_FAILED.value
-
- m0 = fit_summary['g0mean']
- s0 = fit_summary['g0sigma']
- s1 = fit_summary['g1sigma']
- s2 = fit_summary['g2sigma']
- chi2_ndof = fit_summary['chi2_ndof']
- d01 = fit_summary['g1mean'] - m0
-
- mask = 0
- if not fit_summary['is_valid']:
- mask |= BadPixelsFF.FIT_FAILED.value
-
- if not fit_summary['has_accurate_covar']:
- mask |= BadPixelsFF.ACCURATE_COVAR.value
-
- if not is_inside(m0, *peak_lim):
- mask |= BadPixelsFF.NOISE_PEAK_THRESHOLD.value
-
- if not is_inside(d01, *d0_lim):
- mask |= BadPixelsFF.GAIN_THRESHOLD.value
-
- if not is_inside(chi2_ndof, *chi2_lim):
- mask |= BadPixelsFF.CHI2_THRESHOLD.value
-
- if not (is_inside(s1, *(peak_width_lim[0] * s0)) & is_inside(s2, *(
- peak_width_lim[1] * s0))):
- mask |= BadPixelsFF.PEAK_WIDTH_THRESHOLD.value
-
- return mask
-
-
-def mod_name(modno):
- return f"Q{modno // 4 + 1}M{modno % 4 + 1}"
diff --git a/cal_tools/cal_tools/enums.py b/cal_tools/cal_tools/enums.py
index a516dee2f55dbd8ece632c00f523d58a80b2f578..513b25cb80f2827614987a27eab91ab1c9826d59 100644
--- a/cal_tools/cal_tools/enums.py
+++ b/cal_tools/cal_tools/enums.py
@@ -28,21 +28,6 @@ class BadPixels(Enum):
NON_LIN_RESPONSE_REGION = 0b100000000000000000000 # bit 21
-class BadPixelsFF(Enum):
- """ The SLopesFF Bad Pixel Encoding
- """
-
- FIT_FAILED = 0b000000000000000000001 # bit 1
- CHI2_THRESHOLD = 0b000000000000000000010 # bit 2
- NOISE_PEAK_THRESHOLD = 0b000000000000000000100 # bit 3
- GAIN_THRESHOLD = 0b000000000000000001000 # bit 4
- PEAK_WIDTH_THRESHOLD = 0b000000000000000010000 # bit 5
- ACCURATE_COVAR = 0b000000000000000100000 # bit 6
- BAD_DARK = 0b000000000000001000000 # bit 6
- NO_ENTRY = 0b000000000000010000000 # bit 7
- GAIN_DEVIATION = 0b000000000000100000000 # bit 8
-
-
class SnowResolution(Enum):
""" An Enum specifying how to resolve snowy pixels
"""
diff --git a/notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_NBC.ipynb b/notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_NBC.ipynb
index 9ee0f347e3fae0f7219e407d51205cc7ec633d69..6c7a36b8f44b6dfabe902d4f38a5ac7504249481 100644
--- a/notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_NBC.ipynb
+++ b/notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_NBC.ipynb
@@ -4,69 +4,59 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Gain Characterization #\n",
- "\n"
+ "# Gain Characterization (Flat Fields) #\n",
+ "\n",
+ "The following code characterizes the gain of the AGIPD detector from flat field data, i.e. data with X-rays of defined intensity. This data should fullfil the following requirements:\n",
+ "\n",
+ "* intensity should be such that single photon peaks are visible\n",
+ "* data for all modules should be present\n",
+ "* no shadowing should occur on any of the modules\n",
+ "* each pixel should have at minimum arround 100 events per photon peak per memory cell\n",
+ "* if central beam edges are visible, they should not be significantly more intense\n",
+ "\n",
+ "Characterization is done by a weighted average algorithm, which evaluates the peak locations for all pixels\n",
+ "and memory cells of a given module. These locations are then fitted to a linear function of the average peak\n",
+ "location in each module, such that it yield a relative gain correction."
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T07:34:24.467244Z",
+ "start_time": "2019-09-13T07:34:24.457261Z"
+ }
+ },
"outputs": [],
"source": [
- "in_folder = \"/gpfs/exfel/exp/SPB/202030/p900138/scratch/karnem/r0203_r0204_v01/\" # the folder to read histograms from, required\n",
- "out_folder = \"/gpfs/exfel/exp/SPB/202030/p900138/scratch/karnem/r0203_r0204_v01/\" # the folder to output to, required\n",
- "hist_file_template = \"hists_m{:02d}_sum.h5\" # the template to use to access histograms\n",
- "modules = [10] # modules to correct, set to -1 for all, range allowed\n",
- "\n",
- "image_data_path = \"/gpfs/exfel/exp/MID/202030/p900137/raw\" # Path to image data used to create histograms\n",
- "run = 449 # of the run of image data used to create histograms\n",
- "\n",
- "karabo_id = \"MID_DET_AGIPD1M-1\" # karabo karabo_id\n",
- "karabo_da = ['-1'] # a list of data aggregators names, Default [-1] for selecting all data aggregators\n",
- "receiver_id = \"{}CH0\" # inset for receiver devices\n",
- "path_template = 'RAW-R{:04d}-{}-S{:05d}.h5' # the template to use to access data\n",
- "h5path = 'INSTRUMENT/{}/DET/{}:xtdf/' # path in the HDF5 file to images\n",
- "h5path_idx = 'INDEX/{}/DET/{}:xtdf/' # path in the HDF5 file to images\n",
- "h5path_ctrl = '/CONTROL/{}/MDL/FPGA_COMP' # path to control information\n",
- "karabo_id_control = \"MID_IRU_AGIPD1M1\" # karabo-id for control device\n",
- "karabo_da_control = 'AGIPD1MCTRL00' # karabo DA for control infromation\n",
- "\n",
- "use_dir_creation_date = True # use the creation data of the input dir for database queries\n",
- "cal_db_interface = \"tcp://max-exfl016:8015#8045\" # the database interface to use\n",
- "cal_db_timeout = 30000 # in milli seconds\n",
+ "# the following lines should be adjusted depending on data\n",
+ "in_folder = '/gpfs/exfel/exp/MID/201931/p900091/raw/' # path to input data, required\n",
+ "modules = [3] # modules to work on, required, range allowed\n",
+ "out_folder = \"/gpfs/exfel/exp/MID/201931/p900091/usr/FF/2.2\" # path to output to, required\n",
+ "runs = [484, 485,] # runs to use, required, range allowed\n",
+ "sequences = [0,1,2,3]#,4,5,6,7,8] #,5,6,7,8,9,10] # sequences files to use, range allowed\n",
+ "cluster_profile = \"noDB\" # The ipcluster profile to use\n",
"local_output = True # output constants locally\n",
"db_output = False # output constants to database\n",
- "\n",
- "# Fit parameters\n",
- "peak_range = [-30, 30, 35, 70, 95, 135, 145, 220] # where to look for the peaks, [a0, b0, a1, b1, ...] exactly 8 elements\n",
- "peak_width_range = [0, 30, 0, 35, 0, 40, 0, 45] # fit limits on the peak widths, [a0, b0, a1, b1, ...] exactly 8 elements\n",
- "peak_norm_range = [0.0, -1, 0, -1, 0, -1, 0, -1] # \n",
- "\n",
- "# Bad-pixel thresholds (gain evaluation error). Contribute to BadPixel bit \"Gain_Evaluation_Error\"\n",
- "peak_lim = [-30, 30] # Limit of position of noise peak\n",
- "d0_lim = [10, 80] # hard limits for distance between noise and first peak\n",
- "peak_width_lim = [0.9, 1.55, 0.95, 1.65] # hard limits on the peak widths for first and second peak, in units of the noise peak. 4 parameters.\n",
- "chi2_lim = [0, 3.0] # Hard limit on chi2/nDOF value\n",
- "\n",
- "intensity_lim = 15 # Threshold on standard deviation of a histogram in ADU. Contribute to BadPixel bit \"No_Entry\"\n",
- "gain_lim = [0.85, 1.15] # Threshold on gain in relative number. Contribute to BadPixel bit \"Gain_deviation\"\n",
- "\n",
- "cell_range = [1, 3] # range of cell to be considered, [0,0] for all\n",
- "pixel_range = [0, 0, 32, 32] # range of pixels x1,y1,x2,y2 to consider [0,0,512,128] for all\n",
- "max_bins = 0 # Maximum number of bins to consider, 0 for all bins\n",
- "batch_size = [1, 8, 8] # batch size: [cell,x,y]\n",
- "fit_range = [0, 0] # range of a histogram considered for fitting in ADU. Dynamically evaluated in case [0,0]\n",
- "n_peaks_fit = 4 # Number of gaussian peaks to fit including noise peak\n",
- "fix_peaks = False # Fix distance between photon peaks\n",
- "do_minos = False # This is additional feature of minuit to evaluate errors. \n",
- "\n",
- "# Detector conditions\n",
- "max_cells = 0 # number of memory cells used, set to 0 to automatically infer\n",
- "bias_voltage = 300 # Bias voltage\n",
- "acq_rate = 0. # the detector acquisition rate, use 0 to try to auto-determine\n",
- "gain_setting = 0.1 # the gain setting, use 0.1 to try to auto-determine\n",
- "photon_energy = 8.05 # photon energy in keV"
+ "bias_voltage = 300 # detector bias voltage\n",
+ "cal_db_interface = \"tcp://max-exfl016:8026#8035\" # the database interface to use\n",
+ "mem_cells = 0 # number of memory cells used\n",
+ "interlaced = False # assume interlaced data format, for data prior to Dec. 2017\n",
+ "fit_hook = True # fit a hook function to medium gain slope\n",
+ "rawversion = 2 # RAW file format version\n",
+ "instrument = \"MID\"\n",
+ "photon_energy = 9.2 # the photon energy in keV\n",
+ "offset_store = \"\" # for file-baed access\n",
+ "high_res_badpix_3d = False # set this to True if you need high-resolution 3d bad pixel plots. Runtime: ~ 1h\n",
+ "db_input = True # retreive data from calibration database, setting offset-store will overwrite this\n",
+ "deviation_threshold = 0.75 # peaks with an absolute location deviation larger than the medium are are considere bad\n",
+ "acqrate = 0. # acquisition rate\n",
+ "use_dir_creation_date = True\n",
+ "creation_time = \"\" # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC.ms e.g. 2019-07-04 11:02:41.00\n",
+ "gain_setting = 0.1 # gain setting can have value 0 or 1, Default=0.1 for no (None) gain-setting\n",
+ "karabo_da_control = \"AGIPD1MCTRL00\" # karabo DA for control infromation\n",
+ "h5path_ctrl = '/CONTROL/{}/MDL/FPGA_COMP_TEST' # path to control information"
]
},
{
@@ -75,43 +65,147 @@
"metadata": {},
"outputs": [],
"source": [
- "import glob\n",
- "from multiprocessing import Pool\n",
- "import os\n",
- "import traceback\n",
+ "# std library imports\n",
+ "from datetime import datetime\n",
+ "import dateutil\n",
+ "from functools import partial\n",
"import warnings\n",
+ "warnings.filterwarnings('ignore')\n",
+ "import os\n",
"\n",
"import h5py\n",
- "from iminuit import Minuit\n",
- "import matplotlib.pyplot as plt\n",
+ "# numpy and matplot lib specific\n",
"import numpy as np\n",
- "import sharedmem\n",
- "from XFELDetAna.plotting.heatmap import heatmapPlot\n",
- "from XFELDetAna.plotting.simpleplot import simplePlot\n",
- "import XFELDetAna.xfelpyanatools as xana\n",
+ "import matplotlib\n",
+ "matplotlib.use(\"Agg\")\n",
+ "import matplotlib.pyplot as plt\n",
+ "%matplotlib inline\n",
"\n",
- "from cal_tools.ana_tools import save_dict_to_hdf5, get_range\n",
- "from cal_tools.agipdutils_ff import (gauss_fun, get_starting_par, any_in,\n",
- " multi_gauss, set_par_limits, get_mask, fit_n_peaks)\n",
- "from cal_tools.enums import BadPixelsFF\n",
+ "# parallel processing via ipcluster\n",
+ "# make sure a cluster is running with ipcluster start --n=32, give it a while to start\n",
+ "from ipyparallel import Client\n",
"\n",
- "# %load_ext autotime\n",
- "%matplotlib inline\n",
- "warnings.filterwarnings('ignore')"
+ "client = Client(profile=cluster_profile)\n",
+ "view = client[:]\n",
+ "view.use_dill()\n",
+ "\n",
+ "# pyDetLib imports\n",
+ "import XFELDetAna.xfelpycaltools as xcal\n",
+ "import XFELDetAna.xfelpyanatools as xana\n",
+ "from XFELDetAna.util import env\n",
+ "env.iprofile = cluster_profile\n",
+ "profile = cluster_profile\n",
+ "\n",
+ "from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions\n",
+ "from cal_tools.agipdlib import get_num_cells, get_acq_rate, get_gain_setting\n",
+ "from cal_tools.enums import BadPixels\n",
+ "from cal_tools.influx import InfluxLogger\n",
+ "from cal_tools.plotting import show_overview, plot_badpix_3d\n",
+ "from cal_tools.tools import gain_map_files, parse_runs, run_prop_seq_from_path, get_notebook_name, get_dir_creation_date, get_random_db_interface\n"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T07:34:24.830847Z",
+ "start_time": "2019-09-13T07:34:24.745094Z"
+ }
+ },
"outputs": [],
"source": [
- "peak_range = np.reshape(peak_range,(4,2))\n",
- "peak_width_range = np.reshape(peak_width_range,(4,2))\n",
- "peak_width_lim = np.reshape(peak_width_lim,(2,2))\n",
- "peak_norm_range = [None if x == -1 else x for x in peak_norm_range]\n",
- "peak_norm_range = np.reshape(peak_norm_range,(4,2))\n",
- "module = modules[0]"
+ "# usually no need to change these lines\n",
+ "sensor_size = [128, 512]\n",
+ "block_size = [128, 512]\n",
+ "QUADRANTS = 4\n",
+ "MODULES_PER_QUAD = 4\n",
+ "DET_FILE_INSET = \"AGIPD\"\n",
+ "\n",
+ "# Define constant creation time.\n",
+ "if creation_time:\n",
+ " try:\n",
+ " creation_time = datetime.strptime(creation_time, '%Y-%m-%d %H:%M:%S.%f')\n",
+ " except Exception as e:\n",
+ " print(f\"creation_time value error: {e}.\" \n",
+ " \"Use same format as YYYY-MM-DD HR:MN:SC.ms e.g. 2019-07-04 11:02:41.00/n\")\n",
+ " creation_time = None\n",
+ " print(\"Given creation time wont be used.\")\n",
+ "else:\n",
+ " creation_time = None\n",
+ "\n",
+ "if not creation_time and use_dir_creation_date:\n",
+ " ntries = 100\n",
+ " while ntries > 0:\n",
+ " try:\n",
+ " creation_time = get_dir_creation_date(in_folder, runs[0])\n",
+ " break\n",
+ " except OSError:\n",
+ " pass\n",
+ " ntries -= 1\n",
+ " \n",
+ "print(\"Using {} as creation time\".format(creation_time))\n",
+ " \n",
+ "runs = parse_runs(runs)\n",
+ "\n",
+ "if offset_store != \"\":\n",
+ " db_input = False\n",
+ "else:\n",
+ " db_input = True \n",
+ "\n",
+ "limit_trains = 20\n",
+ "limit_trains_eval = None\n",
+ "\n",
+ "print(\"Parameters are:\")\n",
+ "\n",
+ "print(\"Runs: {}\".format(runs))\n",
+ "print(\"Modules: {}\".format(modules))\n",
+ "print(\"Sequences: {}\".format(sequences))\n",
+ "print(\"Using DB: {}\".format(db_output))\n",
+ "\n",
+ "\n",
+ "if instrument == \"SPB\":\n",
+ " loc = \"SPB_DET_AGIPD1M-1\"\n",
+ " dinstance = \"AGIPD1M1\"\n",
+ " karabo_id_control = \"SPB_IRU_AGIPD1M1\"\n",
+ "else:\n",
+ " loc = \"MID_DET_AGIPD1M-1\"\n",
+ " dinstance = \"AGIPD1M2\"\n",
+ " karabo_id_control = \"MID_EXP_AGIPD1M1\"\n",
+ "\n",
+ "cal_db_interface = get_random_db_interface(cal_db_interface)\n",
+ "\n",
+ "# these lines can usually stay as is\n",
+ "fbase = \"{}/{{}}/RAW-{{}}-AGIPD{{:02d}}-S{{:05d}}.h5\".format(in_folder)\n",
+ "gains = np.arange(3)\n",
+ "\n",
+ "\n",
+ "run, prop, seq = run_prop_seq_from_path(in_folder)\n",
+ "logger = InfluxLogger(detector=\"AGIPD\", instrument=instrument, mem_cells=mem_cells,\n",
+ " notebook=get_notebook_name(), proposal=prop)\n",
+ "\n",
+ "# extract the memory cells and acquisition rate from \n",
+ "# the file of the first module, first sequence and first run\n",
+ "channel = 0\n",
+ "fname = fbase.format(runs[0], runs[0].upper(), channel, sequences[0])\n",
+ "if acqrate == 0.:\n",
+ " acqrate = get_acq_rate((fname, loc, channel))\n",
+ " \n",
+ "\n",
+ "if mem_cells == 0:\n",
+ " cells = get_num_cells(fname, loc, channel)\n",
+ " mem_cells = cells # avoid setting twice\n",
+ " \n",
+ "\n",
+ "IL_MODE = interlaced\n",
+ "max_cells = mem_cells if not interlaced else mem_cells*2\n",
+ "memory_cells = mem_cells\n",
+ "print(\"Interlaced mode: {}\".format(IL_MODE))\n",
+ "\n",
+ "cells = np.arange(max_cells)\n",
+ "\n",
+ "print(f\"Acquisition rate and memory cells set from file {fname} are \" \n",
+ " f\"{acqrate} MHz and {max_cells}, respectively\")"
]
},
{
@@ -120,386 +214,912 @@
"metadata": {},
"outputs": [],
"source": [
- "def idx_gen(batch_start, batch_size):\n",
- " \"\"\"\n",
- " This generator iterate across pixels and memory cells starting\n",
- " from batch_start until batch_start+batch_size\n",
- " \"\"\"\n",
- " for c_idx in range(batch_start[0], batch_start[0]+batch_size[0]):\n",
- " for x_idx in range(batch_start[1], batch_start[1]+batch_size[1]):\n",
- " for y_idx in range(batch_start[2], batch_start[2]+batch_size[2]):\n",
- " yield(c_idx, x_idx, y_idx)"
+ "control_fname = f'{in_folder}/{runs[0]}/RAW-{runs[0].upper()}-{karabo_da_control}-S00000.h5'\n",
+ "\n",
+ "if \"{\" in h5path_ctrl:\n",
+ " h5path_ctrl = h5path_ctrl.format(karabo_id_control)\n",
+ "\n",
+ "if gain_setting == 0.1:\n",
+ " if creation_time.replace(tzinfo=None) < dateutil.parser.parse('2020-01-31'):\n",
+ " print(\"Set gain-setting to None for runs taken before 2020-01-31\")\n",
+ " gain_setting = None\n",
+ " else:\n",
+ " try:\n",
+ " gain_setting = get_gain_setting(control_fname, h5path_ctrl)\n",
+ " except Exception as e:\n",
+ " print(f'Error while reading gain setting from: \\n{control_fname}')\n",
+ " print(e)\n",
+ " print(\"Gain setting is not found in the control information\")\n",
+ " print(\"Data will not be processed\")\n",
+ " sequences = []\n",
+ "print(f\"Gain setting: {gain_setting}\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "For the characterization offset maps for each module are needed. The are read in either locally or through the XFEL calibration database."
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T07:34:45.876476Z",
+ "start_time": "2019-09-13T07:34:25.656979Z"
+ }
+ },
"outputs": [],
"source": [
- "n_pixels_x = pixel_range[2]-pixel_range[0]\n",
- "n_pixels_y = pixel_range[3]-pixel_range[1]\n",
+ "from dateutil import parser\n",
+ "offset_g = {}\n",
+ "noise_g = {}\n",
+ "thresholds_g = {}\n",
+ "pc_g = {}\n",
+ "if not db_input:\n",
+ " print(\"Offset, noise and thresholds have been read in from: {}\".format(offset_store))\n",
+ " store_file = h5py.File(offset_store, \"r\")\n",
+ " for i in modules:\n",
+ " qm = \"Q{}M{}\".format(i//4+1, i%4+1)\n",
+ " offset_g[qm] = np.array(store_file[\"{}/Offset/0/data\".format(qm)])\n",
+ " noise_g[qm] = np.array(store_file[\"{}/Noise/0/data\".format(qm)])\n",
+ " thresholds_g[qm] = np.array(store_file[\"{}/Threshold/0/data\".format(qm)])\n",
+ " store_file.close()\n",
+ "else:\n",
+ " print(\"Offset, noise and thresholds have been read in from calibration database:\")\n",
+ " first_ex = True\n",
+ " for i in modules:\n",
+ " qm = \"Q{}M{}\".format(i//4+1, i%4+1)\n",
+ " metadata = ConstantMetaData()\n",
+ " offset = Constants.AGIPD.Offset()\n",
+ " metadata.calibration_constant = offset\n",
+ " det = getattr(Detectors, dinstance)\n",
"\n",
- "hist_data = {}\n",
- "with h5py.File(f\"{in_folder}/{hist_file_template.format(module)}\", 'r') as hf:\n",
- " hist_data['cellId'] = np.array(hf['cellId'][()])\n",
- " hist_data['hRange'] = np.array(hf['hRange'][()])\n",
- " hist_data['nBins'] = np.array(hf['nBins'][()])\n",
- " \n",
- " if cell_range == [0,0]:\n",
- " cell_range[1] = hist_data['cellId'].shape[0]\n",
- " \n",
- " if max_bins == 0:\n",
- " max_bins = hist_data['nBins']\n",
- " \n",
- " hist_data['cellId'] = hist_data['cellId'][cell_range[0]:cell_range[1]]\n",
- " hist_data['hist'] = np.array(hf['hist'][cell_range[0]:cell_range[1], :max_bins, :])\n",
+ " # set the operating condition\n",
+ " condition = Conditions.Dark.AGIPD(memory_cells=max_cells, bias_voltage=bias_voltage,\n",
+ " acquisition_rate=acqrate, gain_setting=gain_setting)\n",
+ " metadata.detector_condition = condition\n",
"\n",
- "n_cells = cell_range[1]-cell_range[0]\n",
- "hist_data['hist'] = hist_data['hist'].reshape(n_cells, max_bins, 512, 128)\n",
- "hist_data['hist'] = hist_data['hist'][:,:,pixel_range[0]:pixel_range[2],pixel_range[1]:pixel_range[3]]\n",
+ " # specify the a version for this constant\n",
+ " if creation_time is None:\n",
+ " metadata.calibration_constant_version = Versions.Now(device=getattr(det, qm))\n",
+ " else:\n",
+ " metadata.calibration_constant_version = Versions.Timespan(device=getattr(det, qm), start=creation_time)\n",
+ " metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=300000)\n",
+ " offset_g[qm] = offset.data\n",
+ " if first_ex:\n",
+ " print(\"Offset map was injected on: {}\".format(metadata.calibration_constant_version.begin_at))\n",
+ " \n",
+ " metadata = ConstantMetaData()\n",
+ " noise = Constants.AGIPD.Noise()\n",
+ " metadata.calibration_constant = noise\n",
"\n",
- "print(f'Data shape {hist_data[\"hist\"].shape}')\n",
- " \n",
- "bin_edges = np.linspace(hist_data['hRange'][0], hist_data['hRange'][1], int(hist_data['nBins']+1))\n",
- "x = (bin_edges[1:] + bin_edges[:-1])[:max_bins] * 0.5\n",
- " \n",
+ " # set the operating condition\n",
+ " condition = Conditions.Dark.AGIPD(memory_cells=max_cells, bias_voltage=bias_voltage,\n",
+ " acquisition_rate=acqrate, gain_setting=gain_setting)\n",
+ " metadata.detector_condition = condition\n",
"\n",
- "work = []\n",
- "for c_idx in range(0, n_cells, batch_size[0]):\n",
- " for x_idx in range(0, n_pixels_x, batch_size[1]):\n",
- " for y_idx in range(0, n_pixels_y, batch_size[2]):\n",
- " work.append([c_idx,x_idx,y_idx])\n",
+ " # specify the a version for this constant\n",
+ " if creation_time is None:\n",
+ " metadata.calibration_constant_version = Versions.Now(device=getattr(det, qm))\n",
+ " else:\n",
+ " metadata.calibration_constant_version = Versions.Timespan(device=getattr(det, qm), start=creation_time)\n",
+ " metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=300000)\n",
+ " noise_g[qm] = noise.data\n",
+ " if first_ex:\n",
+ " print(\"Noise map was injected on: {}\".format(metadata.calibration_constant_version.begin_at))\n",
" \n",
- "print(f'Number of batches {len(work)}')"
+ " metadata = ConstantMetaData()\n",
+ " thresholds = Constants.AGIPD.ThresholdsDark()\n",
+ " metadata.calibration_constant = thresholds\n",
+ "\n",
+ " # set the operating condition\n",
+ " condition = Conditions.Dark.AGIPD(memory_cells=max_cells, bias_voltage=bias_voltage,\n",
+ " acquisition_rate=acqrate, gain_setting=gain_setting)\n",
+ " metadata.detector_condition = condition\n",
+ "\n",
+ " # specify the a version for this constant\n",
+ " if creation_time is None:\n",
+ " metadata.calibration_constant_version = Versions.Now(device=getattr(det, qm))\n",
+ " else:\n",
+ " metadata.calibration_constant_version = Versions.Timespan(device=getattr(det, qm), start=creation_time)\n",
+ " metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=300000)\n",
+ " thresholds_g[qm] = thresholds.data\n",
+ " if first_ex:\n",
+ " print(\"Threshold map was injected on: {}\".format(metadata.calibration_constant_version.begin_at))\n",
+ " \n",
+ " \n",
+ " metadata = ConstantMetaData()\n",
+ " pc = Constants.AGIPD.SlopesPC()\n",
+ " metadata.calibration_constant = pc\n",
+ "\n",
+ " # set the operating condition\n",
+ " condition = Conditions.Dark.AGIPD(memory_cells=max_cells, bias_voltage=bias_voltage,\n",
+ " acquisition_rate=acqrate, gain_setting=gain_setting)\n",
+ " metadata.detector_condition = condition\n",
+ "\n",
+ " # specify the a version for this constant\n",
+ " if creation_time is None:\n",
+ " metadata.calibration_constant_version = Versions.Now(device=getattr(det, qm))\n",
+ " else:\n",
+ " metadata.calibration_constant_version = Versions.Timespan(device=getattr(det, qm), start=creation_time)\n",
+ " metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=300000)\n",
+ " pc_g[qm] = np.nanmedian(pc.data[0,...])/pc.data\n",
+ " if first_ex:\n",
+ " print(\"PC map was injected on: {}\".format(metadata.calibration_constant_version.begin_at))\n",
+ " first_ex = False\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Initial peak estimates ##\n",
+ "\n",
+ "First initial peak estimates need to be defined. This is done by inspecting histograms created from (a subset of) the offset-corrected data for peak locations and peak heights."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
- "scrolled": true
+ "ExecuteTime": {
+ "end_time": "2019-09-13T07:51:37.815384Z",
+ "start_time": "2019-09-13T07:34:45.879121Z"
+ }
},
"outputs": [],
"source": [
- "def fit_banch(batch_start):\n",
- " current_result = {}\n",
- " prev = None\n",
- " for c_idx, x_idx, y_idx in idx_gen(batch_start, batch_size):\n",
- " #print(c_idx, x_idx, y_idx)\n",
+ "def hist_single_module(fbase, runs, sequences, sensor_size, memory_cells, block_size,\n",
+ " il_mode, limit_trains, profile, rawversion, instrument, inp):\n",
+ " \"\"\" This function calculates a per-pixel histogram for a single module\n",
+ " \n",
+ " Runs and sequences give the data to calculate histogram from\n",
+ " \"\"\"\n",
+ " channel, offset, thresholds, pc, noise = inp\n",
+ " \n",
+ " import XFELDetAna.xfelpycaltools as xcal\n",
+ " import numpy as np\n",
+ " import h5py\n",
+ " from XFELDetAna.util import env\n",
+ " \n",
+ " \n",
+ " env.iprofile = profile\n",
+ " \n",
+ " def baseline_correct_via_noise(d, noise, g, baseline_corr_noise_threshold=1000):\n",
+ " \"\"\" Correct baseline shifts by evaluating position of the noise peak\n",
+ "\n",
+ " :param: d: the data to correct, should be a single image\n",
+ " :param noise: the mean noise for the cell id of `d`\n",
+ " :param g: gain array matching `d` array\n",
+ "\n",
+ " Correction is done by identifying the left-most (significant) peak\n",
+ " in the histogram of `d` and shifting it to be centered on 0.\n",
+ " This is done via a continous wavelet transform (CWT), using a Ricker\n",
+ " wavelet.\n",
+ "\n",
+ " Only high gain data is evaulated, and data larger than 50 ADU,\n",
+ " equivalent of roughly a 9 keV photon is ignored.\n",
+ "\n",
+ " Two passes are executed: the first shift is accurate to 10 ADU and\n",
+ " will catch baseline shifts smaller then -2000 ADU. CWT is evaluated\n",
+ " for peaks of widths one, three and five time the noise.\n",
+ " The baseline is then shifted by the position of the summed CWT maximum.\n",
+ "\n",
+ " In a second pass histogram is evaluated within a range of\n",
+ " +- 5*noise of the initial shift value, for peak of width noise.\n",
+ "\n",
+ " Baseline shifts larger than the maximum threshold or positive shifts\n",
+ " larger 10 are discarded, and the original data is returned, otherwise\n",
+ " the shift corrected data is returned.\n",
+ "\n",
+ " \"\"\"\n",
+ " import copy\n",
+ " from scipy.signal import cwt, ricker\n",
+ " # we work on a copy to be able to filter values by setting them to\n",
+ " # np.nan\n",
+ " dd = copy.copy(d)\n",
+ " dd[g != 0] = np.nan # only high gain data\n",
+ " dd[dd > 50] = np.nan # only noise data\n",
+ " h, e = np.histogram(dd, bins=210, range=(-2000, 100))\n",
+ " c = (e[1:] + e[:-1]) / 2\n",
+ "\n",
" try:\n",
- " y = hist_data['hist'][c_idx, :, x_idx, y_idx]\n",
- "\n",
- " if prev is None:\n",
- " prev, _ = get_starting_par(x, y, peak_range, n_peaks=n_peaks_fit)\n",
- "\n",
- " if fit_range == [0, 0]:\n",
- " frange = (prev[f'g0mean']-2*prev[f'g0sigma'],\n",
- " prev[f'g{n_peaks_fit-1}mean'] + prev[f'g{n_peaks_fit-1}sigma'])\n",
- " else:\n",
- " frange = fit_range\n",
- "\n",
- " set_par_limits(prev, peak_range, peak_norm_range,\n",
- " peak_width_range, n_peaks_fit)\n",
- " minuit = fit_n_peaks(x, y, prev, frange,\n",
- " do_minos=do_minos, n_peaks=n_peaks_fit,\n",
- " fix_d01=fix_peaks)\n",
- "\n",
- " ndof = np.rint(frange[1]-frange[0])-len(minuit.args)\n",
- " current_result['chi2_ndof'] = minuit.fval/ndof\n",
- " current_result.update(minuit.fitarg)\n",
- " current_result.update(minuit.get_fmin())\n",
- "\n",
- " fit_result['chi2_ndof'][c_idx, x_idx, y_idx] = current_result['chi2_ndof']\n",
- " fit_result['g0mean'][c_idx, x_idx, y_idx] = current_result['g0mean']\n",
- " fit_result['g1mean'][c_idx, x_idx, y_idx] = current_result['g1mean']\n",
- " fit_result['g2mean'][c_idx, x_idx, y_idx] = current_result['g2mean']\n",
- " fit_result['g3mean'][c_idx, x_idx, y_idx] = current_result['g3mean']\n",
- "\n",
- " for key in minuit.fitarg.keys():\n",
- " if key in fit_result:\n",
- " fit_result[key][c_idx, x_idx, y_idx] = minuit.fitarg[key]\n",
- "\n",
- " fit_result['mask'][c_idx, x_idx, y_idx] = get_mask(current_result,\n",
- " peak_lim,\n",
- " d0_lim, chi2_lim,\n",
- " peak_width_lim)\n",
- " except Exception as e:\n",
- " fit_result['mask'][c_idx, x_idx,\n",
- " y_idx] = BadPixelsFF.FIT_FAILED.value\n",
- " print(c_idx, x_idx, y_idx, e, traceback.format_exc())\n",
+ " cwtmatr = cwt(h, ricker, [noise, 3. * noise, 5. * noise])\n",
+ " except:\n",
+ " return d\n",
+ " cwtall = np.sum(cwtmatr, axis=0)\n",
+ " marg = np.argmax(cwtall)\n",
+ " pc = c[marg]\n",
+ "\n",
+ " high_res_range = (dd > pc - 5 * noise) & (dd < pc + 5 * noise)\n",
+ " dd[~high_res_range] = np.nan\n",
+ " h, e = np.histogram(dd, bins=200,\n",
+ " range=(pc - 5 * noise, pc + 5 * noise))\n",
+ " c = (e[1:] + e[:-1]) / 2\n",
+ " try:\n",
+ " cwtmatr = cwt(h, ricker, [noise, ])\n",
+ " except:\n",
+ " return d\n",
+ " marg = np.argmax(cwtmatr)\n",
+ " pc = c[marg]\n",
+ " # too large shift to be easily decernable via noise\n",
+ " if pc > 10 or pc < -baseline_corr_noise_threshold:\n",
+ " return d\n",
+ " return d - pc\n",
"\n",
- " if fit_result['mask'][c_idx, x_idx, y_idx] == 0:\n",
- " prev = minuit.fitarg\n",
+ " \n",
+ " # function needs to be inline for parallell processing\n",
+ " def read_fun(filename, channel):\n",
+ " \"\"\" A reader function used by pyDetLib\n",
+ " \"\"\"\n",
+ " infile = h5py.File(filename, \"r\", driver=\"core\")\n",
+ " if rawversion == 2:\n",
+ " count = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/count\".format(instrument, channel)])\n",
+ " first = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/first\".format(instrument, channel)])\n",
+ " last_index = int(first[count != 0][-1]+count[count != 0][-1])\n",
+ " first_index = int(first[count != 0][0])\n",
" else:\n",
- " prev = None"
+ " status = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/status\".format(instrument, channel)])\n",
+ " if np.count_nonzero(status != 0) == 0:\n",
+ " return\n",
+ " last = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/last\".format(instrument, channel)])\n",
+ " last_index = int(last[status != 0][-1])\n",
+ " first_index = int(last[status != 0][0])\n",
+ " if limit_trains is not None:\n",
+ " last_index = min(limit_trains*memory_cells+first_index, last_index)\n",
+ " \n",
+ " im = np.array(infile[\"/INSTRUMENT/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/data\".format(instrument, channel)][first_index:last_index,...]) \n",
+ " carr = infile[\"/INSTRUMENT/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/cellId\".format(instrument, channel)][first_index:last_index]\n",
+ " cells = np.squeeze(np.array(carr))\n",
+ " infile.close()\n",
+ " \n",
+ " if il_mode:\n",
+ " ga = im[1::2, 0, ...]\n",
+ " im = im[0::2, 0, ...].astype(np.float32)\n",
+ " else:\n",
+ " ga = im[:, 1, ...]\n",
+ " im = im[:, 0, ...].astype(np.float32)\n",
+ "\n",
+ " im = np.rollaxis(im, 2)\n",
+ " im = np.rollaxis(im, 2, 1)\n",
+ "\n",
+ " ga = np.rollaxis(ga, 2)\n",
+ " ga = np.rollaxis(ga, 2, 1)\n",
+ " return im, ga, cells\n",
+ " \n",
+ " offset_cor = xcal.OffsetCorrection(sensor_size,\n",
+ " offset,\n",
+ " nCells=memory_cells,\n",
+ " blockSize=block_size,\n",
+ " gains=[0,1,2])\n",
+ " offset_cor.mapper = offset_cor._view.map_sync\n",
+ " offset_cor.debug() # force non-parallel processing since outer function will run concurrently\n",
+ " hist_calc = xcal.HistogramCalculator(sensor_size,\n",
+ " bins=4000,\n",
+ " range=(-4000, 8000),\n",
+ " blockSize=block_size)\n",
+ " hist_calc.mapper = hist_calc._view.map_sync\n",
+ " hist_calc.debug() # force non-parallel processing since outer function will run concurrently\n",
+ " for run in runs:\n",
+ " for seq in sequences:\n",
+ " fname = fbase.format(run, run.upper(), channel, seq)\n",
+ " d, ga, c = read_fun(fname, channel)\n",
+ " vidx = (c < offset.shape[2]) & (c < thresholds.shape[2])\n",
+ " c = c[vidx]\n",
+ " d = d[:,:,vidx]\n",
+ " ga = ga[:,:,vidx]\n",
+ " # we need to do proper gain thresholding\n",
+ " g = np.zeros(ga.shape, np.uint8)\n",
+ " g[...] = 2\n",
+ " \n",
+ " g[ga < thresholds[...,c,1]] = 1\n",
+ " g[ga < thresholds[...,c,0]] = 0\n",
+ " d = offset_cor.correct(d, cellTable=c, gainMap=g)\n",
+ " for i in range(d.shape[2]):\n",
+ " mn_noise = np.nanmean(noise[...,c[i]])\n",
+ " d[...,i] = baseline_correct_via_noise(d[...,i],\n",
+ " mn_noise,\n",
+ " g[..., i])\n",
+ " \n",
+ " d *= np.moveaxis(pc[c,...], 0, 2)\n",
+ " \n",
+ " hist_calc.fill(d)\n",
+ " h, e, c, _ = hist_calc.get()\n",
+ " return h, e, c\n",
+ "\n",
+ "inp = []\n",
+ "for i in modules:\n",
+ " qm = \"Q{}M{}\".format(i//4+1, i%4+1)\n",
+ " inp.append((i, offset_g[qm], thresholds_g[qm], pc_g[qm][0,...], noise_g[qm][...,0]))\n",
+ " \n",
+ "p = partial(hist_single_module, fbase, runs, sequences,\n",
+ " sensor_size, memory_cells, block_size, IL_MODE, limit_trains, profile, rawversion, instrument)\n",
+ "res_uncorr = list(map(p, inp))\n"
]
},
{
- "cell_type": "markdown",
- "metadata": {},
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T07:51:38.024035Z",
+ "start_time": "2019-09-13T07:51:37.818276Z"
+ },
+ "scrolled": false
+ },
+ "outputs": [],
"source": [
- "# Single fit ##\n",
+ "d = []\n",
+ "qms = []\n",
+ "for i, r in enumerate(res_uncorr):\n",
+ " ii = list(modules)[i]\n",
+ " qm = \"Q{}M{}\".format(ii//4+1, ii%4+1)\n",
+ " qms.append(qm)\n",
+ " h, e, c = r \n",
+ " d.append({\n",
+ " 'x': c, \n",
+ " 'y': h,\n",
+ " 'drawstyle': 'steps-mid'\n",
+ " })\n",
+ " \n",
+ "fig = xana.simplePlot(d, y_log=False,\n",
+ " figsize=\"2col\",\n",
+ " aspect=2,\n",
+ " x_range=(-50, 500),\n",
+ " x_label=\"Intensity (ADU)\",\n",
+ " y_label=\"Counts\")\n",
"\n",
- "Left plot shows starting parameters for fitting. Right plot shows result of the fit. Errors are evaluated with minos."
+ "fig.savefig(\"{}/FF_module_{}_peak_pos.png\".format(out_folder, \",\".join(qms)))"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T07:51:38.029704Z",
+ "start_time": "2019-09-13T07:51:38.026111Z"
+ }
+ },
"outputs": [],
"source": [
- "hist = hist_data['hist'][1,:,1, 1]\n",
- "prev, shapes = get_starting_par(x, hist, peak_range, n_peaks=n_peaks_fit)\n",
+ "# these should be quite stable\n",
"\n",
- "if fit_range == [0, 0]:\n",
- " frange = (prev[f'g0mean']-2*prev[f'g0sigma'],\n",
- " prev[f'g3mean'] + prev[f'g3sigma'])\n",
- "else:\n",
- " frange = fit_range\n",
- "\n",
- "set_par_limits(prev, peak_range, peak_norm_range,\n",
- " peak_width_range, n_peaks=n_peaks_fit)\n",
- "minuit = fit_n_peaks(x, hist, prev, frange,\n",
- " do_minos=True, n_peaks=n_peaks_fit,\n",
- " fix_d01=fix_peaks)\n",
- "\n",
- "\n",
- "res = minuit.fitarg\n",
- "err = minuit.errors\n",
- "p = minuit.args\n",
- "ya = np.arange(0,1e4)\n",
- "y = multi_gauss(x,n_peaks_fit, *p)\n",
- "peak_colors = ['g', 'y', 'b', 'orange']\n",
- "\n",
- "d = [{'x': x,\n",
- " 'y': hist.astype(np.float64),\n",
- " 'y_err': np.sqrt(hist),\n",
- " 'drawstyle': 'bars',\n",
- " 'errorstyle': 'bars',\n",
- " 'ecolor': 'navy',\n",
- " 'errorcoarsing': 3,\n",
- " 'label': 'X-ray Data'\n",
- " },\n",
- " {'x': x,\n",
- " 'y': y,\n",
- " 'y2': (hist-y)/np.sqrt(hist),\n",
- " 'color': 'red',\n",
- " 'drawstyle':'line',\n",
- " 'drawstyle2': 'steps-mid',\n",
- " 'label': 'Fit'\n",
- " },\n",
- " ]\n",
- "\n",
- "for i in range(n_peaks_fit):\n",
- " d.append({'x': x,\n",
- " 'y': multi_gauss(x, 1, res[f'g{i}n'], res[f'g{i}mean'], res[f'g{i}sigma']),\n",
- " 'drawstyle':'line',\n",
- " 'color': peak_colors[i],\n",
- " })\n",
- " d.append({'x': np.full_like(ya, res[f'g{i}mean']),\n",
- " 'y': ya,\n",
- " 'drawstyle': 'line',\n",
- " 'linestyle': 'dashed',\n",
- " 'color': peak_colors[i],\n",
- " 'label': f'peak {i} = {res[f\"g{i}mean\"]:0.1f} $ \\pm $ {err[f\"g{i}mean\"]:0.2f} ADU' })\n",
- "\n",
- "fig = plt.figure(figsize=(16,7))\n",
- "ax = fig.add_subplot(121)\n",
- "for i, shape in enumerate(shapes):\n",
- " idx = shape[3]\n",
- " plt.errorbar(x[idx], hist[idx], np.sqrt(hist[idx]),\n",
- " marker='+', ls=''\n",
- " )\n",
- " yg = gauss_fun(x[idx], *shape[:3])\n",
- " l = f'Peak {i}: {shape[1]:0.1f} $ \\pm $ {shape[2]:0.2f} ADU'\n",
- " plt.plot(x[idx], yg, label=l)\n",
- "plt.grid(True)\n",
- "plt.xlabel(\"Signal [ADU]\")\n",
- "plt.ylabel(\"Counts\")\n",
- "#plt.ylim(0, np.max(hist)*1.4)\n",
- "plt.legend(ncol=2)\n",
- "\n",
- "\n",
- "ax2 = fig.add_subplot(122)\n",
- "fig2 = xana.simplePlot(d, \n",
- " use_axis=ax2,\n",
- " x_label='Signal [ADU]', \n",
- " y_label='Counts',\n",
- " secondpanel=True, y_log=False, \n",
- " x_range=(frange[0], frange[1]), \n",
- " y_range=(1., np.max(hist)*1.6), \n",
- " legend='top-left-frame-ncol2')\n",
- "\n",
- "print (minuit.get_fmin())\n",
- "minuit.print_matrix()\n",
- "print(minuit.get_param_states())"
+ "peak_estimates = [0, 55, 105, 155]\n",
+ "peak_heights = [5e7, 5e6, 1e6, 5e5]\n",
+ "peak_sigma = [5., 5.,5., 5.,]\n",
+ "\n",
+ "print(\"Using the following peak estimates: {}\".format(peak_estimates))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true
+ },
+ "source": [
+ "## Calculate relative gain per module ##\n",
+ "\n",
+ "Using the so obtained estimates, we calculate the relative gain per module. For this we use the weighted average method implemented in pyDetLib."
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T08:16:50.652492Z",
+ "start_time": "2019-09-13T07:51:38.031787Z"
+ }
+ },
"outputs": [],
"source": [
- "# Allocate memory for fit results\n",
- "fit_result = {}\n",
- "keys = list(minuit.fitarg.keys())\n",
- "keys = [x for x in keys if 'limit_' not in x and 'fix_' not in x]\n",
- "keys += ['chi2_ndof', 'mask', 'gain']\n",
- "for key in keys:\n",
- " dtype = 'f4'\n",
- " if key == 'mask':\n",
- " dtype = 'i4'\n",
- " fit_result[key] = sharedmem.empty([n_cells, n_pixels_x, n_pixels_y], dtype=dtype)"
+ "block_size = [64, 64]\n",
+ "def relgain_single_module(fbase, runs, sequences, peak_estimates,\n",
+ " peak_heights, peak_sigma, memory_cells, sensor_size,\n",
+ " block_size, il_mode, profile, limit_trains_eval, rawversion, instrument, inp):\n",
+ " \"\"\" A function for calculated the relative gain of a single AGIPD module\n",
+ " \"\"\"\n",
+ " \n",
+ " # import needed inline for parallel processing\n",
+ " import XFELDetAna.xfelpycaltools as xcal\n",
+ " import numpy as np\n",
+ " import h5py\n",
+ " from XFELDetAna.util import env\n",
+ " env.iprofile = profile\n",
+ " \n",
+ " channel, offset, thresholds, noise, pc = inp\n",
+ " \n",
+ " def baseline_correct_via_noise(d, noise, g, baseline_corr_noise_threshold=1000):\n",
+ " \"\"\" Correct baseline shifts by evaluating position of the noise peak\n",
+ "\n",
+ " :param: d: the data to correct, should be a single image\n",
+ " :param noise: the mean noise for the cell id of `d`\n",
+ " :param g: gain array matching `d` array\n",
+ "\n",
+ " Correction is done by identifying the left-most (significant) peak\n",
+ " in the histogram of `d` and shifting it to be centered on 0.\n",
+ " This is done via a continous wavelet transform (CWT), using a Ricker\n",
+ " wavelet.\n",
+ "\n",
+ " Only high gain data is evaulated, and data larger than 50 ADU,\n",
+ " equivalent of roughly a 9 keV photon is ignored.\n",
+ "\n",
+ " Two passes are executed: the first shift is accurate to 10 ADU and\n",
+ " will catch baseline shifts smaller then -2000 ADU. CWT is evaluated\n",
+ " for peaks of widths one, three and five time the noise.\n",
+ " The baseline is then shifted by the position of the summed CWT maximum.\n",
+ "\n",
+ " In a second pass histogram is evaluated within a range of\n",
+ " +- 5*noise of the initial shift value, for peak of width noise.\n",
+ "\n",
+ " Baseline shifts larger than the maximum threshold or positive shifts\n",
+ " larger 10 are discarded, and the original data is returned, otherwise\n",
+ " the shift corrected data is returned.\n",
+ "\n",
+ " \"\"\"\n",
+ " import copy\n",
+ " from scipy.signal import cwt, ricker\n",
+ " # we work on a copy to be able to filter values by setting them to\n",
+ " # np.nan\n",
+ " dd = copy.copy(d)\n",
+ " dd[g != 0] = np.nan # only high gain data\n",
+ " dd[dd > 50] = np.nan # only noise data\n",
+ " h, e = np.histogram(dd, bins=210, range=(-2000, 100))\n",
+ " c = (e[1:] + e[:-1]) / 2\n",
+ "\n",
+ " try:\n",
+ " cwtmatr = cwt(h, ricker, [noise, 3. * noise, 5. * noise])\n",
+ " except:\n",
+ " return d\n",
+ " cwtall = np.sum(cwtmatr, axis=0)\n",
+ " marg = np.argmax(cwtall)\n",
+ " pc = c[marg]\n",
+ "\n",
+ " high_res_range = (dd > pc - 5 * noise) & (dd < pc + 5 * noise)\n",
+ " dd[~high_res_range] = np.nan\n",
+ " h, e = np.histogram(dd, bins=200,\n",
+ " range=(pc - 5 * noise, pc + 5 * noise))\n",
+ " c = (e[1:] + e[:-1]) / 2\n",
+ " try:\n",
+ " cwtmatr = cwt(h, ricker, [noise, ])\n",
+ " except:\n",
+ " return d\n",
+ " marg = np.argmax(cwtmatr)\n",
+ " pc = c[marg]\n",
+ " # too large shift to be easily decernable via noise\n",
+ " if pc > 10 or pc < -baseline_corr_noise_threshold:\n",
+ " return d\n",
+ " return d - pc\n",
+ "\n",
+ " \n",
+ " # function needs to be inline for parallell processing\n",
+ " def read_fun(filename, channel):\n",
+ " \"\"\" A reader function used by pyDetLib\n",
+ " \"\"\"\n",
+ " infile = h5py.File(filename, \"r\", driver=\"core\")\n",
+ " if rawversion == 2:\n",
+ " count = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/count\".format(instrument, channel)])\n",
+ " first = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/first\".format(instrument, channel)])\n",
+ " last_index = int(first[count != 0][-1]+count[count != 0][-1])\n",
+ " first_index = int(first[count != 0][0])\n",
+ " else:\n",
+ " status = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/status\".format(instrument, channel)])\n",
+ " if np.count_nonzero(status != 0) == 0:\n",
+ " return\n",
+ " last = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/last\".format(instrument, channel)])\n",
+ " last_index = int(last[status != 0][-1])\n",
+ " first_index = int(last[status != 0][0])\n",
+ " if limit_trains is not None:\n",
+ " last_index = min(limit_trains*memory_cells+first_index, last_index)\n",
+ " im = np.array(infile[\"/INSTRUMENT/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/data\".format(instrument, channel)][first_index:last_index,...]) \n",
+ " carr = infile[\"/INSTRUMENT/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/cellId\".format(instrument, channel)][first_index:last_index]\n",
+ " cells = np.squeeze(np.array(carr))\n",
+ " infile.close()\n",
+ "\n",
+ " if il_mode:\n",
+ " ga = im[1::2, 0, ...]\n",
+ " im = im[0::2, 0, ...].astype(np.float32)\n",
+ " else:\n",
+ " ga = im[:, 1, ...]\n",
+ " im = im[:, 0, ...].astype(np.float32)\n",
+ " \n",
+ " im = np.rollaxis(im, 2)\n",
+ " im = np.rollaxis(im, 2, 1)\n",
+ "\n",
+ " ga = np.rollaxis(ga, 2)\n",
+ " ga = np.rollaxis(ga, 2, 1)\n",
+ " return im, ga, cells\n",
+ " \n",
+ " offset_cor = xcal.OffsetCorrection(sensor_size, offset, nCells=memory_cells,\n",
+ " blockSize=block_size, gains=[0,1,2])\n",
+ " offset_cor.mapper = offset_cor._view.map_sync\n",
+ "\n",
+ " rel_gain = xcal.GainMapCalculator(sensor_size,\n",
+ " peak_estimates,\n",
+ " peak_sigma,\n",
+ " nCells=memory_cells,\n",
+ " peakHeights = peak_heights,\n",
+ " noiseMap=noise,\n",
+ " deviationType=\"relative\")\n",
+ " rel_gain.mapper = rel_gain._view.map_sync\n",
+ " for run in runs:\n",
+ " for seq in sequences:\n",
+ " fname = fbase.format(run, run.upper(), channel, seq)\n",
+ " d, ga, c = read_fun(fname, channel)\n",
+ " # we need to do proper gain thresholding\n",
+ " vidx = (c < offset.shape[2]) & (c < thresholds.shape[2])\n",
+ " c = c[vidx]\n",
+ " d = d[:,:,vidx]\n",
+ " ga = ga[:,:,vidx]\n",
+ " # we need to do proper gain thresholding\n",
+ " g = np.zeros(ga.shape, np.uint8)\n",
+ " g[...] = 2\n",
+ " \n",
+ " g[ga < thresholds[...,c,1]] = 1\n",
+ " g[ga < thresholds[...,c,0]] = 0\n",
+ " d = offset_cor.correct(d, cellTable=c, gainMap=g)\n",
+ " \n",
+ " for i in range(d.shape[2]):\n",
+ " mn_noise = np.nanmean(noise[...,c[i]])\n",
+ " d[...,i] = baseline_correct_via_noise(d[...,i],\n",
+ " mn_noise,\n",
+ " g[..., i])\n",
+ " \n",
+ " d *= np.moveaxis(pc[c,...], 0, 2)\n",
+ " rel_gain.fill(d, cellTable=c)\n",
+ " \n",
+ " gain_map = rel_gain.get()\n",
+ " gain_map_func = rel_gain.getUsingFunc(inverse=False)\n",
+ " \n",
+ " pks, stds = rel_gain.getRawPeaks()\n",
+ " return gain_map, pks, stds, gain_map_func\n",
+ "\n",
+ "inp = []\n",
+ "for i in modules:\n",
+ " qm = \"Q{}M{}\".format(i//4+1, i%4+1)\n",
+ " inp.append((i, offset_g[qm], thresholds_g[qm], noise_g[qm][...,0], pc_g[qm][0,...]))\n",
+ "start = datetime.now()\n",
+ "p = partial(relgain_single_module, fbase, runs, sequences,\n",
+ " peak_estimates, peak_heights, peak_sigma, memory_cells,\n",
+ " sensor_size, block_size, IL_MODE, profile, limit_trains_eval, rawversion, instrument)\n",
+ "res_gain = list(map(p, inp)) # don't run concurently as inner function are parelllized\n",
+ "duration = (datetime.now()-start).total_seconds()\n",
+ "logger.runtime_summary_entry(success=True, runtime=duration)\n",
+ "logger.send()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Finally, we inspect the results, by plotting the number of entries, peak locations and resulting gain maps for each peak. In the course of doing so, we identify bad pixels by either having 0 entries for a peak, or having `nan` values for the peak location."
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T08:16:54.261782Z",
+ "start_time": "2019-09-13T08:16:50.657594Z"
+ }
+ },
"outputs": [],
"source": [
- "# Perform fitting\n",
- "with Pool() as pool:\n",
- " const_out = pool.map(fit_banch, work)"
+ "\n",
+ "gain_m = {}\n",
+ "flatsff = {}\n",
+ "gainoff_g = {}\n",
+ "entries_g = {}\n",
+ "peaks_g = {}\n",
+ "mask_g = {}\n",
+ "cell_to_preview = 4\n",
+ "masks_eval_peaks = [1, 2]\n",
+ "global_eval_peaks = [1]\n",
+ "global_meds = {}\n",
+ "\n",
+ "for i, r in enumerate(res_gain):\n",
+ " ii = list(modules)[i]\n",
+ " qm = \"Q{}M{}\".format(ii//4+1, ii%4+1)\n",
+ " gain, pks, std, gfunc = r\n",
+ " gains, errors, chisq, valid, max_dev, stats = gfunc\n",
+ " _, entries, stds, sow = gain\n",
+ " gain_db = np.zeros((gains.shape[0], gains.shape[1], memory_cells))\n",
+ " gain_mdb = np.zeros((gains.shape[0], gains.shape[1], memory_cells))\n",
+ " entries_db = np.zeros((gains.shape[0], gains.shape[1], memory_cells, 5)) \n",
+ " gainoff_db = np.zeros((gains.shape[0], gains.shape[1], memory_cells))\n",
+ " mask_db = np.zeros((gains.shape[0], gains.shape[1], memory_cells), np.uint32)\n",
+ " \n",
+ " gainoff_g[qm] = gainoff_db\n",
+ " gain_m[qm] = gain_mdb\n",
+ " entries_g[qm] = entries_db\n",
+ " peaks_g[qm] = pks\n",
+ " \n",
+ " # create a mask for unregular pixels\n",
+ " # first bit set if first peak has nan entries\n",
+ " for pk in masks_eval_peaks:\n",
+ " mask_db[~(np.isfinite(gain_mdb) & np.isfinite(gainoff_db))] |= BadPixels.FF_GAIN_EVAL_ERROR.value\n",
+ " mask_db[(np.abs(1-gain_mdb/np.nanmedian(gain_mdb) > deviation_threshold) )] |= BadPixels.FF_GAIN_DEVIATION.value\n",
+ " # second bit set if entries are 0 for first peak\n",
+ " mask_db[entries[...,1] == 0] |= BadPixels.FF_NO_ENTRIES.value\n",
+ " zero_oc = np.count_nonzero((entries > 0).astype(np.int), axis=3)\n",
+ " mask_db[zero_oc <= len(peak_estimates)-2] |= BadPixels.FF_NO_ENTRIES.value \n",
+ " # third bit set if entries of a given adc show significant noise\n",
+ " stds = []\n",
+ " for ii in range(8):\n",
+ " for jj in range(8):\n",
+ " stds.append(np.std(entries_db[ii*16:(ii+1)*16,jj*64+2:(jj+1)*64-2,:,1], axis=(0,1)))\n",
+ " avg_stds = np.median(np.array(stds), axis=0)\n",
+ " \n",
+ " for ii in range(8):\n",
+ " for jj in range(8):\n",
+ " std = np.std(entries_db[ii*16:(ii+1)*16,jj*64+2:(jj+1)*64-2,:,1], axis=(0,1))\n",
+ " if np.any(std > 5*avg_stds):\n",
+ " mask_db[ii*16:(ii+1)*16,jj*64:(jj+1)*64,std > 5*avg_stds] |= BadPixels.FF_GAIN_DEVIATION\n",
+ " \n",
+ " \n",
+ " mask_g[qm] = mask_db\n",
+ " \n",
+ " flat = np.zeros((gains.shape[0], gains.shape[1], memory_cells, 3))\n",
+ " for j in range(2,len(peak_estimates)):\n",
+ " flat[...,j-2] = np.mean(entries[...,j]/np.mean(entries[...,j]))\n",
+ " flat = np.mean(flat, axis=3)\n",
+ " #flat_db = np.zeros((gains.shape[0], gains.shape[1], memory_cells))\n",
+ " #for j in range(2):\n",
+ " # flat_db[...,j::2] = flat\n",
+ " flatsff[qm] = flat\n",
+ " \n",
+ " global_meds[qm] = np.nanmedian(pks[...,global_eval_peaks][np.max(mask_db, axis=2) != 0])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Evaluated peak locations ##\n",
+ "\n",
+ "The following plot shows the evaluated peak locations for each peak. Peak ids increase downward:\n"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T08:16:55.562017Z",
+ "start_time": "2019-09-13T08:16:54.264843Z"
+ },
+ "scrolled": false
+ },
"outputs": [],
"source": [
- "# Evaluate bad pixels\n",
- "fit_result['gain'] = (fit_result['g1mean'] - fit_result['g0mean'])/photon_energy\n",
- "\n",
- "# Calculate histogram width and evaluate cut\n",
- "h_sums = np.sum(hist_data['hist'], axis=1)\n",
- "hist_norm = hist_data['hist'] / h_sums[:, None, :, :]\n",
- "hist_mean = np.sum(hist_norm[:, :max_bins, ...] *\n",
- " x[None, :, None, None], axis=1)\n",
- "hist_sqr = (x[None, :, None, None] - hist_mean[:, None, ...])**2\n",
- "hist_std = np.sqrt(np.sum(hist_norm[:, :max_bins, ...] * hist_sqr, axis=1))\n",
- "\n",
- "fit_result['mask'][hist_std<intensity_lim] |= BadPixelsFF.NO_ENTRY.value\n",
- "\n",
- "# Bad pixel on gain deviation\n",
- "gains = np.copy(fit_result['gain'])\n",
- "gains[fit_result['mask']>0] = np.nan\n",
- "gain_mean = np.nanmean(gains, axis=(1,2))\n",
- "\n",
- "fit_result['mask'][fit_result['gain'] > gain_mean[:,None,None]*gain_lim[1] ] |= BadPixelsFF.GAIN_DEVIATION.value\n",
- "fit_result['mask'][fit_result['gain'] < gain_mean[:,None,None]*gain_lim[0] ] |= BadPixelsFF.GAIN_DEVIATION.value\n"
+ "from mpl_toolkits.axes_grid1 import AxesGrid\n",
+ "cell_to_preview = 4\n",
+ "for qm, data in peaks_g.items():\n",
+ " print(\"The module shown is: {}\".format(qm))\n",
+ " print(\"The cell shown is: {}\".format(cell_to_preview))\n",
+ " print(\"Evaluated peaks: {}\".format(data.shape[-1]))\n",
+ " fig = plt.figure(figsize=(15,15))\n",
+ " grid = AxesGrid(fig, 111, \n",
+ " nrows_ncols=(data.shape[-1], 1),\n",
+ " axes_pad=0.1,\n",
+ " share_all=True,\n",
+ " label_mode=\"L\",\n",
+ " cbar_location=\"right\",\n",
+ " cbar_mode=\"each\",\n",
+ " cbar_size=\"7%\",\n",
+ " cbar_pad=\"2%\")\n",
+ " \n",
+ " for j in range(data.shape[-1]):\n",
+ " d = data[...,cell_to_preview,j]\n",
+ " d[~np.isfinite(d)] = 0\n",
+ " \n",
+ " im = grid[j].imshow(d, interpolation=\"nearest\", vmin=0.9*np.nanmedian(d), vmax=max(1.1*np.nanmedian(d), 50))\n",
+ " cb = grid.cbar_axes[j].colorbar(im)\n",
+ " cb.set_label_text(\"Peak location (ADU)\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Gain Slopes And Offsets ##\n",
+ "\n",
+ "The gain slopes and offsets are deduced by fitting a linear function $y = mx + b$ to the evaluated peaks. Gains are normalized to all pixels and all memory cells of a module by using the average peak locations a $x$ values. Thus slopes centered around 1 are to be expected.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T08:16:56.879582Z",
+ "start_time": "2019-09-13T08:16:55.564572Z"
+ },
+ "scrolled": false
+ },
"outputs": [],
"source": [
- "# Save fit results\n",
- "os.makedirs(out_folder, exist_ok=True)\n",
- "out_name = f'{out_folder}/fits_m{module:02d}.h5'\n",
- "print(f'Save to file: {out_name}')\n",
- "save_dict_to_hdf5({'data': fit_result}, out_name)"
+ "for i, r in enumerate(res_gain):\n",
+ " ii = list(modules)[i]\n",
+ " qm = \"Q{}M{}\".format(ii//4+1, ii%4+1)\n",
+ " gain, pks, std, gfunc = r\n",
+ " gains, errors, chisq, valid, max_dev, stats = gfunc \n",
+ " _, entries, stds, sow = gain\n",
+ " \n",
+ " fig = plt.figure(figsize=(15,8))\n",
+ " ax = fig.add_subplot(211)\n",
+ " im = ax.imshow(gains[...,0], interpolation=\"nearest\", vmin=0.85, vmax=1.15)\n",
+ " cb = fig.colorbar(im)\n",
+ " cb.set_label(\"Gain slope m\")\n",
+ " fig.savefig(\"{}/gain_m_mod{}.png\".format(out_folder, qm))\n",
+ " \n",
+ " ax = fig.add_subplot(212)\n",
+ " ax.hist(gains[...,0].flatten(), range=(0.5, 1.5), bins=100)\n",
+ " ax.set_ylabel(\"Occurences\")\n",
+ " ax.set_xlabel(\"Gain slope m\")\n",
+ " ax.semilogy()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Summary across cells ##"
+ "The gain offsets b are expected to be centered around 0 and minimal, as data is already offset substracted."
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T08:16:58.113910Z",
+ "start_time": "2019-09-13T08:16:56.881383Z"
+ }
+ },
"outputs": [],
"source": [
- "labels = ['Noise peak [ADU]',\n",
- " 'First photon peak [ADU]',\n",
- " f\"gain [ADU/keV], $\\gamma$={photon_energy} [keV]\",\n",
- " \"$\\chi^2$/nDOF\",\n",
- " 'Fraction of bad pixels']\n",
- "for i, key in enumerate(['g0mean', 'g1mean', 'gain', 'chi2_ndof', 'mask']):\n",
- "\n",
- " fig = plt.figure(figsize=(20,5))\n",
- " ax = fig.add_subplot(121)\n",
- " data = fit_result[key]\n",
- " if key == 'mask':\n",
- " data = data>0\n",
- " vmin, vmax = [0, 1]\n",
- " else:\n",
- " vmin, vmax = get_range(data, 5)\n",
- " _ = heatmapPlot(np.mean(data, axis=0).T,\n",
- " add_panels=False, cmap='viridis', use_axis=ax,\n",
- " vmin=vmin, vmax=vmax, lut_label=labels[i] )\n",
- " \n",
- " if key != 'mask':\n",
- " vmin, vmax = get_range(data, 7)\n",
- " ax1 = fig.add_subplot(122)\n",
- " _ = xana.histPlot(ax1,data.flatten(), \n",
- " bins=45,range=[vmin, vmax],\n",
- " log=True,color='red',histtype='stepfilled')\n",
- " plt.xlabel(labels[i])\n",
- " plt.ylabel(\"Counts\") "
+ "for i, r in enumerate(res_gain):\n",
+ " ii = list(modules)[i]\n",
+ " qm = \"Q{}M{}\".format(ii//4+1, ii%4+1)\n",
+ " gain, pks, std, gfunc = r\n",
+ " gains, errors, chisq, valid, max_dev, stats = gfunc\n",
+ " _, entries, stds, sow = gain\n",
+ " \n",
+ " fig = plt.figure(figsize=(15,8))\n",
+ " ax = fig.add_subplot(211)\n",
+ " \n",
+ " im = ax.imshow(gains[...,1], interpolation=\"nearest\", vmin=-25, vmax=25)\n",
+ " cb = fig.colorbar(im)\n",
+ " cb.set_label(\"Gain offset b\")\n",
+ " fig.savefig(\"{}/gain_b_mod{}.png\".format(out_folder, qm))\n",
+ " \n",
+ " ax = fig.add_subplot(212)\n",
+ " ax.hist(gains[...,1].flatten(), range=(-25, 25), bins=100)\n",
+ " ax.set_ylabel(\"Occurences\")\n",
+ " ax.set_xlabel(\"Gain offset b\")\n",
+ " ax.semilogy()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## histograms of fit parameters ##"
+ "## Bad Pixels ##\n",
+ "\n",
+ "Bad pixels are defined as any of the following criteria: \n",
+ "\n",
+ "* the gain evaluation did not converge, or location of noise peak deviates by more than than the deviation threshold from the median location.\n",
+ "* the number of entries for the noise peak in 0.\n",
+ "* the standard deviation of the number of entries is larger than 5 times the standard deviation for the ASIC the pixel is on."
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T08:18:08.277881Z",
+ "start_time": "2019-09-13T08:18:08.272390Z"
+ }
+ },
"outputs": [],
"source": [
- "fig = plt.figure(figsize=(10, 5))\n",
- "ax0 = fig.add_subplot(111)\n",
- "a = ax0.hist(hist_std.flatten(), bins=100, range=(0,100) )\n",
- "ax0.plot([intensity_lim, intensity_lim], [0, np.nanmax(a[0])], linewidth=1.5, color='red' ) \n",
- "ax0.set_xlabel('Histogram width [ADU]', fontsize=14)\n",
- "ax0.set_ylabel('Number of histograms', fontsize=14)\n",
- "ax0.set_title(f'{hist_std[hist_std<intensity_lim].shape[0]} histograms below threshold in {intensity_lim} ADU',\n",
- " fontsize=14, fontweight='bold')\n",
- "ax0.grid()\n",
- "plt.yscale('log')"
+ "print(\"The deviation threshold is: {}\".format(deviation_threshold))"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T08:18:10.840520Z",
+ "start_time": "2019-09-13T08:18:09.769451Z"
+ }
+ },
"outputs": [],
"source": [
- "num_bins = frange[1] - frange[0]\n",
- "fig = plt.figure(figsize=(16,5))\n",
- "ax = fig.add_subplot(131)\n",
- "_ = xana.histPlot(ax,fit_result['g1mean'].flatten(), \n",
- " bins= num_bins,range=[frange[0] ,frange[1]],\n",
- " log=True,color='red',histtype='stepfilled') \n",
- "plt.xlabel(\"g1 mean [ADU]\")\n",
- "\n",
- "ax1 = fig.add_subplot(132)\n",
- "_ = xana.histPlot(ax1,fit_result['g2mean'].flatten(), \n",
- "# bins=45,range=[80 ,140],\n",
- " bins= num_bins,range=[frange[0] ,frange[1]],\n",
- " log=True,color='red',histtype='stepfilled') \n",
- "plt.xlabel(\"g2 mean [ADU]\")\n",
- "\n",
- "ax2 = fig.add_subplot(133)\n",
- "_ = xana.histPlot(ax2,fit_result['g3mean'].flatten(), \n",
- " bins= num_bins,range=[frange[0] ,frange[1]],\n",
- " log=True,color='red',histtype='stepfilled') \n",
- "_ = plt.xlabel(\"g3 mean [ADU]\")"
+ "for i, r in enumerate(res_gain):\n",
+ " ii = list(modules)[i]\n",
+ " qm = \"Q{}M{}\".format(ii//4+1, ii%4+1)\n",
+ " mask_db = mask_g[qm]\n",
+ " \n",
+ " fig = plt.figure(figsize=(15,8))\n",
+ " ax = fig.add_subplot(211)\n",
+ " im = ax.imshow(np.log2(mask_db[...,cell_to_preview]), interpolation=\"nearest\", vmin=0, vmax=32)\n",
+ " cb = fig.colorbar(im)\n",
+ " cb.set_label(\"Mask value\")\n",
+ " fig.savefig(\"{}/mask_mod{}.png\".format(out_folder, qm))\n",
+ " \n",
+ " ax = fig.add_subplot(212)\n",
+ " ax.hist(np.log2(mask_db.flatten()), range=(0, 32), bins=32, normed=True)\n",
+ " ax.set_ylabel(\"Occurences\")\n",
+ " ax.set_xlabel(\"Mask value\")\n",
+ " ax.semilogy()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T08:18:14.349244Z",
+ "start_time": "2019-09-13T08:18:11.849053Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "cols = {BadPixels.FF_GAIN_EVAL_ERROR.value: (BadPixels.FF_GAIN_EVAL_ERROR.name, '#FF000080'),\n",
+ " BadPixels.FF_GAIN_DEVIATION.value: (BadPixels.FF_GAIN_DEVIATION.name, '#0000FF80'),\n",
+ " BadPixels.FF_NO_ENTRIES.value: (BadPixels.FF_NO_ENTRIES.name, '#00FF0080'),\n",
+ " BadPixels.FF_GAIN_EVAL_ERROR.value |\n",
+ " BadPixels.FF_GAIN_DEVIATION.value: ('EVAL+DEV', '#DD00DD80'),\n",
+ " BadPixels.FF_GAIN_EVAL_ERROR.value |\n",
+ " BadPixels.FF_NO_ENTRIES.value: ('EVAL+NO_ENTRY', '#00DDDD80'),\n",
+ " BadPixels.FF_GAIN_DEVIATION.value |\n",
+ " BadPixels.FF_NO_ENTRIES.value: ('DEV+NO_ENTRY', '#DDDD0080'),\n",
+ " }\n",
+ "\n",
+ "rebin = 32 if not high_res_badpix_3d else 2\n",
+ "\n",
+ "gain = 0\n",
+ "for mod, data in mask_g.items():\n",
+ " plot_badpix_3d(data, cols, title=mod, rebin_fac=rebin)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-07-23T18:02:58.599579Z",
+ "start_time": "2019-07-23T18:02:58.316311Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "if local_output:\n",
+ " ofile = \"{}/agipd_gain_store_{}_modules_{}.h5\".format(out_folder,\n",
+ " \"_\".join([str(r) for r in runs]),\n",
+ " \"_\".join([str(m) for m in modules]))\n",
+ " store_file = h5py.File(ofile, \"w\")\n",
+ " for i, r in enumerate(res_gain):\n",
+ " ii = list(modules)[i]\n",
+ " qm = \"Q{}M{}\".format(ii//4+1, ii%4+1)\n",
+ " gain, pks, std, gfunc = r\n",
+ " gains, errors, chisq, valid, max_dev, stats = gfunc\n",
+ " gainmap, entires, stds, sow = gain\n",
+ " store_file[\"/{}/Gain/0/data\".format(qm)] = gains[...,0]\n",
+ " store_file[\"/{}/GainOffset/0/data\".format(qm)] = gains[...,1]\n",
+ " store_file[\"/{}/Flat/0/data\".format(qm)] = flatsff[qm]\n",
+ " store_file[\"/{}/Entries/0/data\".format(qm)] = entires\n",
+ " store_file[\"/{}/BadPixels/0/data\".format(qm)] = mask_g[qm] \n",
+ " store_file.close()"
]
},
{
@@ -508,96 +1128,446 @@
"metadata": {},
"outputs": [],
"source": [
- "fig = plt.figure(figsize=(16,5))\n",
- "ax = fig.add_subplot(131)\n",
- "_ = xana.histPlot(ax,fit_result['g0sigma'].flatten(), \n",
- " bins=45,range=[0 ,50],\n",
- " log=True,color='red',histtype='stepfilled') \n",
- "plt.xlabel(\"g1 sigma [ADU]\")\n",
- "\n",
- "ax1 = fig.add_subplot(132)\n",
- "_ = xana.histPlot(ax1,fit_result['g1sigma'].flatten(), \n",
- " bins=45,range=[0 ,50],\n",
- " log=True,color='red',histtype='stepfilled') \n",
- "plt.xlabel(\"g2 sigma [ADU]\")\n",
- "\n",
- "ax2 = fig.add_subplot(133)\n",
- "_ = xana.histPlot(ax2,fit_result['g2sigma'].flatten(), \n",
- " bins=45,range=[0 ,50],\n",
- " log=True,color='red',histtype='stepfilled') \n",
- "_ = plt.xlabel(\"g3 sigma [ADU]\")"
+ "proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]\n",
+ "file_loc = proposal + ' ' + ' '.join(list(map(str,runs)))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-12T14:49:28.355922Z",
+ "start_time": "2019-09-12T14:49:28.340063Z"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "if db_output:\n",
+ " for i, r in enumerate(res_gain):\n",
+ " ii = list(modules)[i]\n",
+ " qm = \"Q{}M{}\".format(ii//4+1, ii%4+1)\n",
+ " \n",
+ " gain, pks, std, gfunc = r\n",
+ " gains, errors, chisq, valid, max_dev, stats = gfunc\n",
+ " gainmap, entires, stds, sow = gain\n",
+ " \n",
+ " det = getattr(Detectors, dinstance)\n",
+ " device = getattr(det, qm)\n",
+ " # gains related\n",
+ " metadata = ConstantMetaData()\n",
+ " cgain = Constants.AGIPD.SlopesFF()\n",
+ " cgain.data = gains\n",
+ " metadata.calibration_constant = cgain\n",
+ "\n",
+ " # set the operating condition\n",
+ " condition = Conditions.Illuminated.AGIPD(memory_cells, bias_voltage, 9.2,\n",
+ " pixels_x=512, pixels_y=128, beam_energy=None,\n",
+ " acquisition_rate=acqrate, gain_setting=gain_setting)\n",
+ " \n",
+ " metadata.detector_condition = condition\n",
+ "\n",
+ " # specify the a version for this constant\n",
+ " if creation_time is None:\n",
+ " metadata.calibration_constant_version = Versions.Now(device=getattr(det, qm))\n",
+ " else:\n",
+ " metadata.calibration_constant_version = Versions.Timespan(device=getattr(det, qm), start=creation_time)\n",
+ "\n",
+ " metadata.calibration_constant_version.raw_data_location = file_loc\n",
+ " metadata.send(cal_db_interface, timeout=300000)\n",
+ " \n",
+ " \n",
+ " # bad pixels \n",
+ " metadata = ConstantMetaData()\n",
+ " bp = Constants.AGIPD.BadPixelsFF()\n",
+ " bp.data = mask_g[qm] \n",
+ " metadata.calibration_constant = bp\n",
+ "\n",
+ " # set the operating condition\n",
+ " condition = Conditions.Illuminated.AGIPD(memory_cells, bias_voltage, 9.2,\n",
+ " pixels_x=512, pixels_y=128, beam_energy=None,\n",
+ " acquisition_rate=acqrate, gain_setting=gain_setting)\n",
+ " \n",
+ " metadata.detector_condition = condition\n",
+ "\n",
+ " # specify the a version for this constant\n",
+ " if creation_time is None:\n",
+ " metadata.calibration_constant_version = Versions.Now(device=getattr(det, qm))\n",
+ " else:\n",
+ " metadata.calibration_constant_version = Versions.Timespan(device=getattr(det, qm), start=creation_time)\n",
+ " metadata.calibration_constant_version.raw_data_location = file_loc\n",
+ " metadata.send(cal_db_interface, timeout=300000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Summary across pixels ##\n",
+ "## Sanity check ##\n",
"\n",
- "Mean and median values are calculated across all pixels for each memory cell. "
+ "Finally, we perform a correction of the data used to derive the gain constants with said constants. We expected an improvement in peak FWHM, if constants have been deduced correctly."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T09:23:54.999086Z",
+ "start_time": "2019-09-13T09:16:57.293565Z"
+ },
"scrolled": false
},
"outputs": [],
"source": [
- "def plot_error_band(key, x, ax):\n",
- " \n",
- " cdata = np.copy(fit_result[key])\n",
- " cdata[fit_result['mask']>0] = np.nan\n",
- " \n",
- " mean = np.nanmean(cdata, axis=(1,2))\n",
- " median = np.nanmedian(cdata, axis=(1,2))\n",
- " std = np.nanstd(cdata, axis=(1,2))\n",
- " mad = np.nanmedian(np.abs(cdata - median[:,None,None]), axis=(1,2))\n",
- "\n",
- " ax0 = fig.add_subplot(111)\n",
- " ax0.plot(x, mean, 'k', color='#3F7F4C', label=\" mean value \")\n",
- " ax0.plot(x, median, 'o', color='red', label=\" median value \")\n",
- " ax0.fill_between(x, mean-std, mean+std,\n",
- " alpha=0.6, edgecolor='#3F7F4C', facecolor='#7EFF99',\n",
- " linewidth=1, linestyle='dashdot', antialiased=True,\n",
- " label=\" mean value $ \\pm $ std \")\n",
- "\n",
- " ax0.fill_between(x, median-mad, median+mad,\n",
- " alpha=0.3, edgecolor='red', facecolor='red',\n",
- " linewidth=1, linestyle='dashdot', antialiased=True,\n",
- " label=\" median value $ \\pm $ mad \")\n",
- " \n",
- " if f'error_{key}' in fit_result:\n",
- " cerr = np.copy(fit_result[f'error_{key}'])\n",
- " cerr[fit_result['mask']>0] = np.nan\n",
+ "def hist_single_module(fbase, runs, sequences, il_mode, profile, limit_trains, memory_cells, rawversion, instrument, inp):\n",
+ " channel, offset, thresholds, relgain, noise, pc = inp\n",
+ " gain, pks, std, gfunc = relgain\n",
+ " gains, errors, chisq, valid, max_dev, stats = gfunc\n",
+ " \n",
+ " import XFELDetAna.xfelpycaltools as xcal\n",
+ " import numpy as np\n",
+ " import h5py\n",
+ " import copy\n",
+ " from XFELDetAna.util import env\n",
+ " env.iprofile = profile\n",
+ " sensor_size = [128, 512]\n",
+ " \n",
+ " block_size = sensor_size\n",
+ " \n",
+ " \n",
+ " def baseline_correct_via_noise(d, noise, g, baseline_corr_noise_threshold=1000):\n",
+ " \"\"\" Correct baseline shifts by evaluating position of the noise peak\n",
+ "\n",
+ " :param: d: the data to correct, should be a single image\n",
+ " :param noise: the mean noise for the cell id of `d`\n",
+ " :param g: gain array matching `d` array\n",
+ "\n",
+ " Correction is done by identifying the left-most (significant) peak\n",
+ " in the histogram of `d` and shifting it to be centered on 0.\n",
+ " This is done via a continous wavelet transform (CWT), using a Ricker\n",
+ " wavelet.\n",
+ "\n",
+ " Only high gain data is evaulated, and data larger than 50 ADU,\n",
+ " equivalent of roughly a 9 keV photon is ignored.\n",
+ "\n",
+ " Two passes are executed: the first shift is accurate to 10 ADU and\n",
+ " will catch baseline shifts smaller then -2000 ADU. CWT is evaluated\n",
+ " for peaks of widths one, three and five time the noise.\n",
+ " The baseline is then shifted by the position of the summed CWT maximum.\n",
+ "\n",
+ " In a second pass histogram is evaluated within a range of\n",
+ " +- 5*noise of the initial shift value, for peak of width noise.\n",
+ "\n",
+ " Baseline shifts larger than the maximum threshold or positive shifts\n",
+ " larger 10 are discarded, and the original data is returned, otherwise\n",
+ " the shift corrected data is returned.\n",
+ "\n",
+ " \"\"\"\n",
+ " import copy\n",
+ " from scipy.signal import cwt, ricker\n",
+ " # we work on a copy to be able to filter values by setting them to\n",
+ " # np.nan\n",
+ " dd = copy.copy(d)\n",
+ " dd[g != 0] = np.nan # only high gain data\n",
+ " dd[dd > 50] = np.nan # only noise data\n",
+ " h, e = np.histogram(dd, bins=210, range=(-2000, 100))\n",
+ " c = (e[1:] + e[:-1]) / 2\n",
+ "\n",
+ " try:\n",
+ " cwtmatr = cwt(h, ricker, [noise, 3. * noise, 5. * noise])\n",
+ " except:\n",
+ " return d\n",
+ " cwtall = np.sum(cwtmatr, axis=0)\n",
+ " marg = np.argmax(cwtall)\n",
+ " pc = c[marg]\n",
+ "\n",
+ " high_res_range = (dd > pc - 5 * noise) & (dd < pc + 5 * noise)\n",
+ " dd[~high_res_range] = np.nan\n",
+ " h, e = np.histogram(dd, bins=200,\n",
+ " range=(pc - 5 * noise, pc + 5 * noise))\n",
+ " c = (e[1:] + e[:-1]) / 2\n",
+ " try:\n",
+ " cwtmatr = cwt(h, ricker, [noise, ])\n",
+ " except:\n",
+ " return d\n",
+ " marg = np.argmax(cwtmatr)\n",
+ " pc = c[marg]\n",
+ " # too large shift to be easily decernable via noise\n",
+ " if pc > 10 or pc < -baseline_corr_noise_threshold:\n",
+ " return d\n",
+ " return d - pc\n",
+ " \n",
+ " def read_fun(filename, channel):\n",
+ " \n",
+ " infile = h5py.File(filename, \"r\", driver=\"core\")\n",
+ " if rawversion == 2:\n",
+ " count = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/count\".format(instrument, channel)])\n",
+ " first = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/first\".format(instrument, channel)])\n",
+ " last_index = int(first[count != 0][-1]+count[count != 0][-1])\n",
+ " first_index = int(first[count != 0][0])\n",
+ " else:\n",
+ " status = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/status\".format(instrument, channel)])\n",
+ " if np.count_nonzero(status != 0) == 0:\n",
+ " return\n",
+ " last = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/last\".format(instrument, channel)])\n",
+ " last_index = int(last[status != 0][-1])\n",
+ " first_index = int(last[status != 0][0])\n",
+ " \n",
+ " if limit_trains is not None:\n",
+ " last_index = min(limit_trains*memory_cells+first_index, last_index)\n",
+ " \n",
+ " im = np.array(infile[\"/INSTRUMENT/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/data\".format(instrument, channel)][first_index:last_index,...]) \n",
+ " carr = infile[\"/INSTRUMENT/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/cellId\".format(instrument, channel)][first_index:last_index]\n",
+ " cells = np.squeeze(np.array(carr))\n",
+ " infile.close()\n",
" \n",
- " meanerr = np.nanmean(cerr, axis=(1,2))\n",
- " ax0.fill_between(x, mean-meanerr, mean+meanerr,\n",
- " alpha=0.6, edgecolor='#089FFF', facecolor='#089FFF',\n",
- " linewidth=1, linestyle='dashdot', antialiased=True,\n",
- " label=\" mean fit error \")\n",
" \n",
+ " if il_mode:\n",
+ " ga = im[1::2, 0, ...]\n",
+ " im = im[0::2, 0, ...].astype(np.float32)\n",
+ " else:\n",
+ " ga = im[:, 1, ...]\n",
+ " im = im[:, 0, ...].astype(np.float32)\n",
"\n",
- "x = np.linspace(*cell_range, n_cells)\n",
+ " im = np.rollaxis(im, 2)\n",
+ " im = np.rollaxis(im, 2, 1)\n",
"\n",
- "for i, key in enumerate(['g0mean', 'g1mean', 'gain', 'chi2_ndof']):\n",
+ " ga = np.rollaxis(ga, 2)\n",
+ " ga = np.rollaxis(ga, 2, 1)\n",
+ " return im, ga, cells\n",
+ " \n",
+ " offset_cor = xcal.OffsetCorrection(sensor_size, offset, nCells=memory_cells, blockSize=block_size, gains=[0,1,2])\n",
+ " offset_cor.debug()\n",
+ " \n",
+ " hist_calc = xcal.HistogramCalculator(sensor_size, bins=2000, range=(0, 2000), blockSize=block_size)\n",
+ " hist_calc.debug()\n",
+ "\n",
+ " hist_calc_uncorr = xcal.HistogramCalculator(sensor_size, bins=2000, range=(0, 2000), blockSize=block_size)\n",
+ " hist_calc_uncorr.debug()\n",
+ " \n",
+ " \n",
+ " for run in runs:\n",
+ " for seq in sequences[:2]:\n",
+ " \n",
+ " fname = fbase.format(run, run.upper(), channel, seq)\n",
+ " \n",
+ " d, ga, c = read_fun(fname, channel)\n",
+ " \n",
+ " # we need to do proper gain thresholding\n",
+ " vidx = (c < offset.shape[2]) & (c < thresholds.shape[2])\n",
+ " c = c[vidx]\n",
+ " d = d[:,:,vidx]\n",
+ " ga = ga[:,:,vidx]\n",
"\n",
- " fig = plt.figure(figsize=(10, 5))\n",
- " ax0 = fig.add_subplot(111)\n",
- " plot_error_band(key, x, ax0)\n",
+ " g = np.zeros(ga.shape, np.uint8)\n",
+ " g[...] = 2\n",
+ " \n",
+ " g[ga < thresholds[...,c,1]] = 1\n",
+ " g[ga < thresholds[...,c,0]] = 0\n",
+ " d = offset_cor.correct(d, cellTable=c, gainMap=g)\n",
+ " for i in range(d.shape[2]):\n",
+ " mn_noise = np.nanmean(noise[...,c[i]])\n",
+ " d[...,i] = baseline_correct_via_noise(d[...,i],\n",
+ " mn_noise,\n",
+ " g[..., i])\n",
+ " \n",
+ " d *= np.moveaxis(pc[c,...], 0, 2)\n",
+ " \n",
+ " hist_calc_uncorr.fill(d)\n",
+ " d = (d-gains[..., 1][...,None])/gains[..., 0][...,None]\n",
+ " #d = d/gains[..., 0][...,None]\n",
+ " hist_calc.fill(d) \n",
+ " \n",
+ " h, e, c, _ = hist_calc.get()\n",
+ " hu = hist_calc_uncorr.get()\n",
+ " return h, e, c, hu[0]\n",
"\n",
- " ax0.set_xlabel('Memory Cell ID', fontsize=14)\n",
- " ax0.set_ylabel(labels[i], fontsize=14)\n",
- " ax0.grid()\n",
- " _ = ax0.legend()\n"
+ "inp = []\n",
+ "idx = 0\n",
+ "for i in modules:\n",
+ " qm = \"Q{}M{}\".format(i//4+1, i%4+1)\n",
+ " \n",
+ " inp.append((i, offset_g[qm], thresholds_g[qm], res_gain[idx], noise_g[qm][...,0], pc_g[qm][0,...]))\n",
+ " idx += 1\n",
+ " \n",
+ "p = partial(hist_single_module, fbase, runs, sequences, IL_MODE, profile, limit_trains, memory_cells, rawversion, instrument)\n",
+ "#res = view.map_sync(p, inp)\n",
+ "res = list(map(p, inp))\n"
]
},
{
- "cell_type": "markdown",
- "metadata": {},
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "ExecuteTime": {
+ "end_time": "2019-09-13T09:23:56.023378Z",
+ "start_time": "2019-09-13T09:23:55.002226Z"
+ }
+ },
+ "outputs": [],
"source": [
- "## Cut flow ##"
+ "from iminuit import Minuit\n",
+ "from iminuit.util import make_func_code, describe\n",
+ "from IPython.display import HTML, display\n",
+ "import tabulate\n",
+ "\n",
+ "# fitting\n",
+ "par_ests = {}\n",
+ "par_ests[\"mu0\"] = 0\n",
+ "par_ests[\"mu1\"] = 50\n",
+ "par_ests[\"mu2\"] = 100\n",
+ "par_ests[\"limit_mu0\"] = [-5, 5]\n",
+ "par_ests[\"limit_mu1\"] = [35, 65]\n",
+ "par_ests[\"limit_mu2\"] = [100, 150]\n",
+ "par_ests[\"s0\"] = 10\n",
+ "par_ests[\"s1\"] = 10\n",
+ "par_ests[\"s2\"] = 10\n",
+ "par_ests[\"limit_A0\"] = [1e5, 1e12]\n",
+ "par_ests[\"limit_A1\"] = [1e4, 1e10]\n",
+ "par_ests[\"limit_A2\"] = [1e3, 1e10]\n",
+ "\n",
+ "\n",
+ "par_ests[\"throw_nan\"] = False\n",
+ "par_ests[\"pedantic\"] = False\n",
+ "par_ests[\"print_level\"] = 1\n",
+ "\n",
+ "def gaussian3(x, mu0, s0, A0, mu1, s1, A1, mu2, s2, A2):\n",
+ " return (A0/np.sqrt(2*np.pi*s0**2)*np.exp(-0.5*((x-mu0)/s0)**2) +\n",
+ " A1/np.sqrt(2*np.pi*s1**2)*np.exp(-0.5*((x-mu1)/s1)**2) +\n",
+ " A2/np.sqrt(2*np.pi*s2**2)*np.exp(-0.5*((x-mu2)/s2)**2))\n",
+ "\n",
+ "\n",
+ "f_sig = describe(gaussian3)[1:]\n",
+ "\n",
+ "class _Chi2Functor:\n",
+ " def __init__(self, f, x, y, err):\n",
+ " self.f = f\n",
+ " self.x = x\n",
+ " self.y = y\n",
+ " self.err = err\n",
+ " f_sig = describe(f)\n",
+ " # this is how you fake function\n",
+ " # signature dynamically\n",
+ " self.func_code = make_func_code(\n",
+ " f_sig[1:]) # docking off independent variable\n",
+ " self.func_defaults = None # this keeps numpy.vectorize happy\n",
+ "\n",
+ " def __call__(self, *arg):\n",
+ " # notice that it accept variable length\n",
+ " # positional arguments\n",
+ " # chi2 = sum((y-self.f(x,*arg))**2 for x,y in zip(self.x,self.y))\n",
+ " return np.sum(((self.f(self.x, *arg) - self.y) ** 2) / self.err)\n",
+ " \n",
+ " \n",
+ "d = []\n",
+ "y_range_max = 0\n",
+ "table = []\n",
+ "headers = ['Module',\n",
+ " 'FWHM (cor.) [ADU]', 'Separation (cor.) [$\\sigma$]',\n",
+ " 'FWHM (uncor.) [ADU]', 'Separation (uncor.) [$\\sigma$]', \n",
+ " 'Improvement'\n",
+ " ]\n",
+ "for i, r in enumerate(res):\n",
+ " qm = \"Q{}M{}\".format(i//4+1, i%4+1)\n",
+ " row = []\n",
+ " row.append(qm)\n",
+ " \n",
+ " h, e, c, hu = r\n",
+ " \n",
+ " \n",
+ " d.append({\n",
+ " 'x': c, \n",
+ " 'y': h,\n",
+ " 'drawstyle': 'steps-mid',\n",
+ " 'label': '{}: corrected'.format(qm),\n",
+ " 'marker': '.',\n",
+ " 'color': 'blue',\n",
+ " })\n",
+ " \n",
+ " idx = (h > 0) & np.isfinite(h)\n",
+ " h_non_zero = h[idx]\n",
+ " c_non_zero = c[idx]\n",
+ " par_ests[\"A0\"] = np.float(h[np.argmin(abs(c-0))])\n",
+ " par_ests[\"A1\"] = np.float(h[np.argmin(abs(c-50))])\n",
+ " par_ests[\"A2\"] = np.float(h[np.argmin(abs(c-100))])\n",
+ " wrapped = _Chi2Functor(gaussian3, c_non_zero, h_non_zero,\n",
+ " np.sqrt(h_non_zero))\n",
+ " \n",
+ " m = Minuit(wrapped, **par_ests)\n",
+ " fmin = m.migrad()\n",
+ " \n",
+ " xt = np.arange(0, 200)\n",
+ " \n",
+ " yt = gaussian3(xt, m.values['mu0'], m.values['s0'], m.values['A0'],\n",
+ " m.values['mu1'], m.values['s1'], m.values['A1'],\n",
+ " m.values['mu2'], m.values['s2'], m.values['A2'])\n",
+ "\n",
+ " d.append({\n",
+ " 'x': xt, \n",
+ " 'y': yt,\n",
+ " 'label': '{}: corrected (fit)'.format(qm),\n",
+ " 'color': 'green',\n",
+ " 'drawstyle': 'line',\n",
+ " 'linewidth': 2\n",
+ " })\n",
+ " \n",
+ " \n",
+ " d.append({\n",
+ " 'x': c, \n",
+ " 'y': hu,\n",
+ " 'drawstyle': 'steps-mid',\n",
+ " 'label': '{}: uncorrected'.format(qm),\n",
+ " 'marker': '.',\n",
+ " 'color': 'grey',\n",
+ " 'alpha': 0.5\n",
+ " })\n",
+ " \n",
+ " row += [m.values['s1']*2.35, (m.values['mu1']-m.values['mu0'])/m.values['s1']]\n",
+ " \n",
+ " \n",
+ " idx = (hu > 0) & np.isfinite(hu)\n",
+ " h_non_zero = hu[idx]\n",
+ " c_non_zero = c[idx]\n",
+ " wrapped = _Chi2Functor(gaussian3, c_non_zero, h_non_zero,\n",
+ " np.sqrt(h_non_zero))\n",
+ " \n",
+ " m = Minuit(wrapped, **par_ests)\n",
+ " fmin = m.migrad()\n",
+ " \n",
+ " xt = np.arange(0, 200)\n",
+ " \n",
+ " yt = gaussian3(xt, m.values['mu0'], m.values['s0'], m.values['A0'],\n",
+ " m.values['mu1'], m.values['s1'], m.values['A1'],\n",
+ " m.values['mu2'], m.values['s2'], m.values['A2'])\n",
+ "\n",
+ " d.append({\n",
+ " 'x': xt, \n",
+ " 'y': yt,\n",
+ " 'label': '{}: uncorrected (fit)'.format(qm),\n",
+ " 'color': 'red',\n",
+ " 'linewidth': 2\n",
+ " })\n",
+ " \n",
+ " row += [m.values['s1']*2.35, (m.values['mu1']-m.values['mu0'])/m.values['s1']]\n",
+ " \n",
+ " row.append(\"{:0.2f} %\".format(100*(row[3]/row[1]-1)))\n",
+ " \n",
+ " y_range_max = max(y_range_max, np.max(h[c>25])*1.5)\n",
+ " \n",
+ " # output table\n",
+ " table.append(row)\n",
+ " \n",
+ "fig = xana.simplePlot(d, y_log=False, figsize=\"2col\",\n",
+ " aspect=2,\n",
+ " x_range=(0, 200),\n",
+ " legend='top-right-frame',\n",
+ " y_range=(0, y_range_max),\n",
+ " x_label='Energy (ADU)',\n",
+ " y_label='Counts')\n",
+ "\n",
+ "display(HTML(tabulate.tabulate(table, tablefmt='html', headers=headers,\n",
+ " numalign=\"right\", floatfmt=\"0.2f\")))"
]
},
{
@@ -605,67 +1575,7 @@
"execution_count": null,
"metadata": {},
"outputs": [],
- "source": [
- "fig = plt.figure(figsize=(10, 5))\n",
- "ax = fig.add_subplot(111)\n",
- "\n",
- "n_bars = 8\n",
- "x = np.arange(n_bars)\n",
- "width = 0.3\n",
- "\n",
- "msk = fit_result['mask']\n",
- "n_fits = np.prod(msk.shape)\n",
- "y = [any_in(msk, BadPixelsFF.FIT_FAILED.value),\n",
- " any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value),\n",
- " any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |\n",
- " BadPixelsFF.CHI2_THRESHOLD.value),\n",
- " any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |\n",
- " BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value),\n",
- " any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |\n",
- " BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value |\n",
- " BadPixelsFF.NOISE_PEAK_THRESHOLD.value),\n",
- " any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |\n",
- " BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value |\n",
- " BadPixelsFF.NOISE_PEAK_THRESHOLD.value | BadPixelsFF.PEAK_WIDTH_THRESHOLD.value),\n",
- " any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |\n",
- " BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value |\n",
- " BadPixelsFF.NOISE_PEAK_THRESHOLD.value | BadPixelsFF.PEAK_WIDTH_THRESHOLD.value\n",
- " | BadPixelsFF.NO_ENTRY.value),\n",
- " any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |\n",
- " BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value |\n",
- " BadPixelsFF.NOISE_PEAK_THRESHOLD.value | BadPixelsFF.PEAK_WIDTH_THRESHOLD.value\n",
- " | BadPixelsFF.NO_ENTRY.value| BadPixelsFF.GAIN_DEVIATION.value)\n",
- " ]\n",
- "\n",
- "y2 = [any_in(msk, BadPixelsFF.FIT_FAILED.value),\n",
- " any_in(msk, BadPixelsFF.ACCURATE_COVAR.value),\n",
- " any_in(msk, BadPixelsFF.CHI2_THRESHOLD.value),\n",
- " any_in(msk, BadPixelsFF.GAIN_THRESHOLD.value),\n",
- " any_in(msk, BadPixelsFF.NOISE_PEAK_THRESHOLD.value),\n",
- " any_in(msk, BadPixelsFF.PEAK_WIDTH_THRESHOLD.value),\n",
- " any_in(msk, BadPixelsFF.NO_ENTRY.value),\n",
- " any_in(msk, BadPixelsFF.GAIN_DEVIATION.value)\n",
- " ]\n",
- "\n",
- "y = (1 - np.sum(y, axis=(1,2,3))/n_fits)*100\n",
- "y2 = (1 - np.sum(y2, axis=(1,2,3))/n_fits)*100\n",
- "\n",
- "labels = ['Fit failes',\n",
- " 'Accurate covar',\n",
- " 'Chi2/nDOF',\n",
- " 'Gain',\n",
- " 'Noise peak',\n",
- " 'Peak width',\n",
- " 'No Entry',\n",
- " 'Gain deviation']\n",
- "\n",
- "ax.bar(x, y2, width, label='Only this cut')\n",
- "ax.bar(x, y, width, label='Cut flow')\n",
- "plt.xticks(x, labels, rotation=90)\n",
- "plt.ylim(y[5]-0.5,100)\n",
- "plt.grid(True)\n",
- "plt.legend()"
- ]
+ "source": []
}
],
"metadata": {
diff --git a/notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_Summary.ipynb b/notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_Summary.ipynb
deleted file mode 100644
index 6134769f935f5dde096812d28ca2c8ee982d77d0..0000000000000000000000000000000000000000
--- a/notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_Summary.ipynb
+++ /dev/null
@@ -1,417 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Gain Characterization Summary #\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "in_folder = \"/gpfs/exfel/exp/SPB/202030/p900138/scratch/karnem/r0203_r0204_v01/\" # the folder to read histograms from, required\n",
- "out_folder = \"/gpfs/exfel/exp/SPB/202030/p900138/scratch/karnem/r0203_r0204_v02\" # the folder to output to, required\n",
- "hist_file_template = \"hists_m{:02d}_sum.h5\"\n",
- "modules = [10] # modules to correct, set to -1 for all, range allowed\n",
- "\n",
- "image_data_path = \"/gpfs/exfel/exp/MID/202030/p900137/raw\"\n",
- "run = 449 # runs of image data used to create histograms\n",
- "\n",
- "karabo_id = \"MID_DET_AGIPD1M-1\" # karabo karabo_id\n",
- "karabo_da = ['-1'] # a list of data aggregators names, Default [-1] for selecting all data aggregators\n",
- "receiver_id = \"{}CH0\" # inset for receiver devices\n",
- "path_template = 'RAW-R{:04d}-{}-S{:05d}.h5' # the template to use to access data\n",
- "h5path = 'INSTRUMENT/{}/DET/{}:xtdf/' # path in the HDF5 file to images\n",
- "h5path_idx = 'INDEX/{}/DET/{}:xtdf/' # path in the HDF5 file to images\n",
- "h5path_ctrl = '/CONTROL/{}/MDL/FPGA_COMP' # path to control information\n",
- "karabo_id_control = \"MID_IRU_AGIPD1M1\" # karabo-id for control device\n",
- "karabo_da_control = 'AGIPD1MCTRL00' # karabo DA for control infromation\n",
- "\n",
- "use_dir_creation_date = True # use the creation data of the input dir for database queries\n",
- "cal_db_interface = \"tcp://max-exfl016:8015#8045\" # the database interface to use\n",
- "cal_db_timeout = 30000 # in milli seconds\n",
- "local_output = True # output constants locally\n",
- "db_output = True # output constants to database\n",
- "\n",
- "# Fit parameters\n",
- "peak_range = [-30,30,35,65,80,130,145,200] # where to look for the peaks, [a0, b0, a1, b1, ...] exactly 8 elements\n",
- "peak_width_range = [0, 30, 0, 35, 0, 40, 0, 45] # fit limits on the peak widths, [a0, b0, a1, b1, ...] exactly 8 elements\n",
- "\n",
- "# Bad-pixel thresholds\n",
- "d0_lim = [10, 70] # hard limits for d0 value (distance between noise and first peak)\n",
- "peak_width_lim = [0.97, 1.43, 1.03, 1.57] # hard limits on the peak widths, [a0, b0, a1, b1, ...] in units of the noise peak. 4 parameters.\n",
- "chi2_lim = [0,3.0] # Hard limit on chi2/nDOF value\n",
- "\n",
- "cell_range = [1,5] # range of cell to be considered, [0,0] for all\n",
- "pixel_range = [0,0,512,128] # range of pixels x1,y1,x2,y2 to consider [0,0,512,128] for all\n",
- "max_bins = 250 # Maximum number of bins to consider\n",
- "batch_size = [1,8,8] # batch size: [cell,x,y]\n",
- "n_peaks_fit = 4 # Number of gaussian peaks to fit including noise peak\n",
- "fix_peaks = True # Fix distance between photon peaks\n",
- "\n",
- "\n",
- "# Detector conditions\n",
- "max_cells = 0 # number of memory cells used, set to 0 to automatically infer\n",
- "bias_voltage = 300 # Bias voltage\n",
- "acq_rate = 0. # the detector acquisition rate, use 0 to try to auto-determine\n",
- "gain_setting = 0.1 # the gain setting, use 0.1 to try to auto-determine\n",
- "photon_energy = 9.2 # photon energy in keV"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "import glob\n",
- "from multiprocessing import Pool\n",
- "import re\n",
- "import os\n",
- "import traceback\n",
- "import warnings\n",
- "\n",
- "from dateutil import parser\n",
- "from extra_geom import AGIPD_1MGeometry\n",
- "import h5py\n",
- "from iCalibrationDB import Detectors, Conditions, Constants\n",
- "from IPython.display import HTML, display, Markdown, Latex\n",
- "import matplotlib.pyplot as plt\n",
- "import numpy as np\n",
- "import tabulate\n",
- "from XFELDetAna.plotting.heatmap import heatmapPlot\n",
- "from XFELDetAna.plotting.simpleplot import simplePlot\n",
- "\n",
- "from cal_tools.ana_tools import save_dict_to_hdf5\n",
- "from cal_tools.ana_tools import get_range\n",
- "from cal_tools.agipdlib import (get_num_cells, get_acq_rate, get_gain_setting)\n",
- "from cal_tools.agipdutils_ff import (any_in, mod_name)\n",
- "from cal_tools.tools import get_dir_creation_date, send_to_db\n",
- "from cal_tools.enums import BadPixels, BadPixelsFF\n",
- "\n",
- "#%load_ext autotime\n",
- "%matplotlib inline\n",
- "warnings.filterwarnings('ignore')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Get operation conditions\n",
- "\n",
- "filename = glob.glob(f\"{image_data_path}/r{run:04d}/*-AGIPD[0-1][0-9]-*\")[0]\n",
- "channel = int(re.findall(r\".*-AGIPD([0-9]+)-.*\", filename)[0])\n",
- "control_fname = f'{image_data_path}/r{run:04d}/RAW-R{run:04d}-{karabo_da_control}-S00000.h5'\n",
- "\n",
- "# Evaluate number of memory cells\n",
- "mem_cells = get_num_cells(filename, karabo_id, channel)\n",
- "if mem_cells is None:\n",
- " raise ValueError(f\"No raw images found in {filename}\")\n",
- "\n",
- "# Evaluate aquisition rate\n",
- "acq_rate = get_acq_rate(filename, karabo_id, channel)\n",
- "\n",
- "# Evaluate creation time\n",
- "creation_time = None\n",
- "if use_dir_creation_date:\n",
- " creation_time = get_dir_creation_date(image_data_path, run)\n",
- " \n",
- "# Evaluate gain setting\n",
- "if gain_setting == 0.1:\n",
- " if creation_time.replace(tzinfo=None) < parser.parse('2020-01-31'):\n",
- " print(\"Set gain-setting to None for runs taken before 2020-01-31\")\n",
- " gain_setting = None\n",
- " else:\n",
- " try:\n",
- " gain_setting = get_gain_setting(control_fname, h5path_ctrl)\n",
- " except Exception as e:\n",
- " print(f'Error while reading gain setting from: \\n{control_fname}')\n",
- " print(e)\n",
- " print(\"Set gain settion to 0\")\n",
- " gain_setting = 0\n",
- " \n",
- "# Evaluate detector instance for mapping\n",
- "instrument = karabo_id.split(\"_\")[0]\n",
- "if instrument == \"SPB\":\n",
- " dinstance = \"AGIPD1M1\"\n",
- "else:\n",
- " dinstance = \"AGIPD1M2\"\n",
- " \n",
- "print(f\"Using {creation_time} as creation time\")\n",
- "print(f\"Operating conditions are:\\n• Bias voltage: {bias_voltage}\\n• Memory cells: {mem_cells}\\n\"\n",
- " f\"• Acquisition rate: {acq_rate}\\n• Gain setting: {gain_setting}\\n• Photon Energy: {photon_energy}\\n\")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Load constants for all modules\n",
- "keys = ['g0mean', 'g1mean', 'gain', 'chi2_ndof', 'mask']\n",
- "\n",
- "fit_data = {}\n",
- "labels = {'g0mean': 'Noise peak position [ADU]',\n",
- " 'g1mean': 'First photon peak [ADU]',\n",
- " 'gain': f\"Gain [ADU/keV], $\\gamma$={photon_energy} [keV]\",\n",
- " 'chi2_ndof': '$\\chi^2$/nDOF',\n",
- " 'mask': 'Fraction of bad pixels over cells' }\n",
- "\n",
- "modules = []\n",
- "for mod in range(16):\n",
- " fit_data[mod] = {}\n",
- " try:\n",
- " hf = h5py.File(f'{out_folder}/fits_m{mod:02d}.h5', 'r')\n",
- " shape = hf['data/g0mean'].shape\n",
- " for key in keys:\n",
- " fit_data[mod][key] = hf[f'data/{key}'][()]\n",
- " #print(shape)\n",
- " print(f\"{in_folder}/{hist_file_template.format(mod)}\") \n",
- " modules.append(mod)\n",
- " except Exception as e:\n",
- " err = f\"Error: {e}\\nError traceback: {traceback.format_exc()}\"\n",
- " #print(err)\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Calculate SlopesFF and BadPixels to be send to DB\n",
- "bpmask = {}\n",
- "slopesFF = {}\n",
- "\n",
- "for mod in modules:\n",
- " bpmask[mod] = np.zeros(fit_data[mod]['mask'].shape).astype(np.int32)\n",
- " bpmask[mod][ any_in(fit_data[mod]['mask'], BadPixelsFF.NO_ENTRY.value) ] = BadPixels.FF_NO_ENTRIES.value\n",
- " bpmask[mod][ any_in(fit_data[mod]['mask'], \n",
- " BadPixelsFF.GAIN_DEVIATION.value) ] |= BadPixels.FF_GAIN_DEVIATION.value\n",
- " bpmask[mod][ any_in(fit_data[mod]['mask'], \n",
- " BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |\n",
- " BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value |\n",
- " BadPixelsFF.NOISE_PEAK_THRESHOLD.value | BadPixelsFF.PEAK_WIDTH_THRESHOLD.value) ] |= BadPixels.FF_GAIN_EVAL_ERROR.value\n",
- " \n",
- " # Set value for bad pixel to average across pixels for a given module\n",
- " slopesFF[mod] = np.copy(fit_data[mod]['gain'])\n",
- " slopesFF[mod][fit_data[mod]['mask']>0] = np.nan\n",
- " gain_mean = np.nanmean(slopesFF[mod], axis=(1,2))\n",
- " \n",
- " for i in range(slopesFF[mod].shape[0]):\n",
- " slopesFF[mod][i][ fit_data[mod]['mask'][i]>0 ] = gain_mean[i]\n",
- " "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Send constants to DB\n",
- "def send_const(mod):\n",
- " try:\n",
- " device = getattr(getattr(Detectors, dinstance), mod_name(mod))\n",
- "\n",
- " # gain\n",
- " constant = Constants.AGIPD.SlopesFF()\n",
- " constant.data = np.moveaxis(np.moveaxis(slopesFF[mod],0,2),0,1)\n",
- " send_to_db(device, constant, condition, file_loc,\n",
- " cal_db_interface, creation_time,\n",
- " verbosity=1, timeout=cal_db_timeout)\n",
- "\n",
- " # bad pixels \n",
- " constant_bp = Constants.AGIPD.BadPixelsFF()\n",
- " constant_bp.data = np.moveaxis(np.moveaxis(bpmask[mod],0,2),0,1)\n",
- " send_to_db(device, constant_bp, condition, file_loc,\n",
- " cal_db_interface, creation_time,\n",
- " verbosity=1, timeout=cal_db_timeout)\n",
- " \n",
- " except Exception as e:\n",
- " err = f\"Error: {e}\\nError traceback: {traceback.format_exc()}\"\n",
- " when = None\n",
- "\n",
- "# set the operating condition\n",
- "condition = Conditions.Illuminated.AGIPD(mem_cells, bias_voltage, 9.2,\n",
- " pixels_x=512, pixels_y=128, beam_energy=None,\n",
- " acquisition_rate=acq_rate, gain_setting=gain_setting)\n",
- "\n",
- "# Check, if we have a shape we expect\n",
- "if db_output:\n",
- " if slopesFF[modules[0]].shape == (mem_cells,512,128):\n",
- " with Pool(processes=16) as pool:\n",
- " const_out = pool.map(send_const, modules)\n",
- " else:\n",
- " print(f\"Constants are not sent to the DB because of the shape mismatsh\")\n",
- " print(f\"Expected {(mem_cells,512,128)}, observed {slopesFF[modules[0]].shape}\")\n",
- " \n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# DEfine AGIPD geometry\n",
- "geom = AGIPD_1MGeometry.from_quad_positions(quad_pos=[\n",
- " (-525, 625),\n",
- " (-550, -10),\n",
- " (520, -160),\n",
- " (542.5, 475),\n",
- "])"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "if cell_range==[0,0]:\n",
- " cell_range[1] = shape[0]\n",
- " \n",
- "const_data = {}\n",
- "for key in keys:\n",
- " const_data[key] = np.full((16,shape[0],512,128), np.nan)\n",
- " for i in range(16):\n",
- " if key in fit_data[i]:\n",
- " const_data[key][i,:,pixel_range[0]:pixel_range[2],\n",
- " pixel_range[1]:pixel_range[3]] = fit_data[i][key]\n",
- " \n",
- "const_data['slopesFF'] = np.full((16,shape[0],512,128), np.nan)\n",
- "labels['slopesFF'] = f'slopesFF [ADU/keV], $\\gamma$={photon_energy} [keV]'\n",
- "for i in range(16):\n",
- " if i in slopesFF:\n",
- " const_data['slopesFF'][i,:,pixel_range[0]:pixel_range[2],\n",
- " pixel_range[1]:pixel_range[3]] = slopesFF[i]\n"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Summary across pixels ##"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "scrolled": true
- },
- "outputs": [],
- "source": [
- "for key in const_data.keys():\n",
- " fig = plt.figure(figsize=(20,20))\n",
- " ax = fig.add_subplot(111)\n",
- " if key=='mask':\n",
- " data = np.nanmean(const_data[key]>0, axis=1)\n",
- " vmin, vmax = (0,1)\n",
- " else:\n",
- " data = np.nanmean(const_data[key], axis=1)\n",
- " vmin, vmax = get_range(data, 5)\n",
- " ax = geom.plot_data_fast(data, ax=ax, cmap=\"jet\", vmin=vmin, vmax=vmax, figsize=(20,20))\n",
- " _ = ax.set_title(labels[key])"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Summary across cells ##\n",
- "\n",
- "Good pixels only."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "for key in const_data.keys():\n",
- " data = np.copy(const_data[key])\n",
- " if key=='mask':\n",
- " data = data>0\n",
- " else:\n",
- " data[const_data['mask']>0] = np.nan\n",
- " \n",
- " d = []\n",
- " for i in range(16):\n",
- " d.append({'x': np.arange(data[i].shape[0]),\n",
- " 'y': np.nanmean(data[i], axis=(1,2)),\n",
- " 'drawstyle': 'steps-pre',\n",
- " 'label': f'{i}',\n",
- " 'linewidth': 2,\n",
- " 'linestyle': '--' if i>7 else '-'\n",
- " })\n",
- "\n",
- " fig = plt.figure(figsize=(15, 6))\n",
- " ax = fig.add_subplot(111)\n",
- "\n",
- " _ = simplePlot(d, xrange=(-12, 510),\n",
- " x_label='Memory Cell ID',\n",
- " y_label=labels[key],\n",
- " use_axis=ax,\n",
- " legend='top-left-frame-ncol8',)\n",
- " ylim = ax.get_ylim()\n",
- " ax.set_ylim(ylim[0], ylim[1] + np.abs(ylim[1]-ylim[0])*0.2)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Summary table ##"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "table = []\n",
- "for i in modules:\n",
- " table.append((i,\n",
- " f\"{np.nanmean(slopesFF[i]):0.1f} +- {np.nanstd(slopesFF[i]):0.2f}\",\n",
- " f\"{np.nanmean(bpmask[i]>0)*100:0.1f} ({np.nansum(bpmask[i]>0)})\"\n",
- " ))\n",
- "md = display(Latex(tabulate.tabulate(table, tablefmt='latex',\n",
- " headers=[\"Module\", \"Gain [ADU/keV]\", \"Bad pixels [%(Count)]\"])))"
- ]
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Python 3",
- "language": "python",
- "name": "python3"
- },
- "language_info": {
- "codemirror_mode": {
- "name": "ipython",
- "version": 3
- },
- "file_extension": ".py",
- "mimetype": "text/x-python",
- "name": "python",
- "nbconvert_exporter": "python",
- "pygments_lexer": "ipython3",
- "version": "3.6.7"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 1
-}
diff --git a/notebooks/AGIPD/playground/Characterize_AGIPD_Gain_FlatFields_NBC.ipynb b/notebooks/AGIPD/playground/Characterize_AGIPD_Gain_FlatFields_NBC.ipynb
deleted file mode 100644
index 9f4a3351b56b1857bd6f79465aa94b5249a03bd5..0000000000000000000000000000000000000000
--- a/notebooks/AGIPD/playground/Characterize_AGIPD_Gain_FlatFields_NBC.ipynb
+++ /dev/null
@@ -1,1602 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Gain Characterization (Flat Fields) #\n",
- "\n",
- "The following code characterizes the gain of the AGIPD detector from flat field data, i.e. data with X-rays of defined intensity. This data should fullfil the following requirements:\n",
- "\n",
- "* intensity should be such that single photon peaks are visible\n",
- "* data for all modules should be present\n",
- "* no shadowing should occur on any of the modules\n",
- "* each pixel should have at minimum arround 100 events per photon peak per memory cell\n",
- "* if central beam edges are visible, they should not be significantly more intense\n",
- "\n",
- "Characterization is done by a weighted average algorithm, which evaluates the peak locations for all pixels\n",
- "and memory cells of a given module. These locations are then fitted to a linear function of the average peak\n",
- "location in each module, such that it yield a relative gain correction."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T07:34:24.467244Z",
- "start_time": "2019-09-13T07:34:24.457261Z"
- }
- },
- "outputs": [],
- "source": [
- "# the following lines should be adjusted depending on data\n",
- "in_folder = '/gpfs/exfel/exp/MID/201931/p900091/raw/' # path to input data, required\n",
- "modules = [3] # modules to work on, required, range allowed\n",
- "out_folder = \"/gpfs/exfel/exp/MID/201931/p900091/usr/FF/2.2\" # path to output to, required\n",
- "runs = [484, 485,] # runs to use, required, range allowed\n",
- "sequences = [0,1,2,3]#,4,5,6,7,8] #,5,6,7,8,9,10] # sequences files to use, range allowed\n",
- "cluster_profile = \"noDB\" # The ipcluster profile to use\n",
- "local_output = True # output constants locally\n",
- "db_output = False # output constants to database\n",
- "bias_voltage = 300 # detector bias voltage\n",
- "cal_db_interface = \"tcp://max-exfl016:8026#8035\" # the database interface to use\n",
- "mem_cells = 0 # number of memory cells used\n",
- "interlaced = False # assume interlaced data format, for data prior to Dec. 2017\n",
- "fit_hook = True # fit a hook function to medium gain slope\n",
- "rawversion = 2 # RAW file format version\n",
- "instrument = \"MID\"\n",
- "photon_energy = 9.2 # the photon energy in keV\n",
- "offset_store = \"\" # for file-baed access\n",
- "high_res_badpix_3d = False # set this to True if you need high-resolution 3d bad pixel plots. Runtime: ~ 1h\n",
- "db_input = True # retreive data from calibration database, setting offset-store will overwrite this\n",
- "deviation_threshold = 0.75 # peaks with an absolute location deviation larger than the medium are are considere bad\n",
- "acqrate = 0. # acquisition rate\n",
- "use_dir_creation_date = True\n",
- "creation_time = \"\" # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC.ms e.g. 2019-07-04 11:02:41.00\n",
- "gain_setting = 0.1 # gain setting can have value 0 or 1, Default=0.1 for no (None) gain-setting\n",
- "karabo_da_control = \"AGIPD1MCTRL00\" # karabo DA for control infromation\n",
- "h5path_ctrl = '/CONTROL/{}/MDL/FPGA_COMP_TEST' # path to control information"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# std library imports\n",
- "from datetime import datetime\n",
- "import dateutil\n",
- "from functools import partial\n",
- "import warnings\n",
- "warnings.filterwarnings('ignore')\n",
- "import os\n",
- "\n",
- "import h5py\n",
- "# numpy and matplot lib specific\n",
- "import numpy as np\n",
- "import matplotlib\n",
- "matplotlib.use(\"Agg\")\n",
- "import matplotlib.pyplot as plt\n",
- "%matplotlib inline\n",
- "\n",
- "# parallel processing via ipcluster\n",
- "# make sure a cluster is running with ipcluster start --n=32, give it a while to start\n",
- "from ipyparallel import Client\n",
- "\n",
- "client = Client(profile=cluster_profile)\n",
- "view = client[:]\n",
- "view.use_dill()\n",
- "\n",
- "# pyDetLib imports\n",
- "import XFELDetAna.xfelpycaltools as xcal\n",
- "import XFELDetAna.xfelpyanatools as xana\n",
- "from XFELDetAna.util import env\n",
- "env.iprofile = cluster_profile\n",
- "profile = cluster_profile\n",
- "\n",
- "from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions\n",
- "from cal_tools.agipdlib import get_num_cells, get_acq_rate, get_gain_setting\n",
- "from cal_tools.enums import BadPixels\n",
- "from cal_tools.influx import InfluxLogger\n",
- "from cal_tools.plotting import show_overview, plot_badpix_3d\n",
- "from cal_tools.tools import gain_map_files, parse_runs, run_prop_seq_from_path, get_notebook_name, get_dir_creation_date, get_random_db_interface\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T07:34:24.830847Z",
- "start_time": "2019-09-13T07:34:24.745094Z"
- }
- },
- "outputs": [],
- "source": [
- "# usually no need to change these lines\n",
- "sensor_size = [128, 512]\n",
- "block_size = [128, 512]\n",
- "QUADRANTS = 4\n",
- "MODULES_PER_QUAD = 4\n",
- "DET_FILE_INSET = \"AGIPD\"\n",
- "\n",
- "# Define constant creation time.\n",
- "if creation_time:\n",
- " try:\n",
- " creation_time = datetime.strptime(creation_time, '%Y-%m-%d %H:%M:%S.%f')\n",
- " except Exception as e:\n",
- " print(f\"creation_time value error: {e}.\" \n",
- " \"Use same format as YYYY-MM-DD HR:MN:SC.ms e.g. 2019-07-04 11:02:41.00/n\")\n",
- " creation_time = None\n",
- " print(\"Given creation time wont be used.\")\n",
- "else:\n",
- " creation_time = None\n",
- "\n",
- "if not creation_time and use_dir_creation_date:\n",
- " ntries = 100\n",
- " while ntries > 0:\n",
- " try:\n",
- " creation_time = get_dir_creation_date(in_folder, runs[0])\n",
- " break\n",
- " except OSError:\n",
- " pass\n",
- " ntries -= 1\n",
- " \n",
- "print(\"Using {} as creation time\".format(creation_time))\n",
- " \n",
- "runs = parse_runs(runs)\n",
- "\n",
- "if offset_store != \"\":\n",
- " db_input = False\n",
- "else:\n",
- " db_input = True \n",
- "\n",
- "limit_trains = 20\n",
- "limit_trains_eval = None\n",
- "\n",
- "print(\"Parameters are:\")\n",
- "\n",
- "print(\"Runs: {}\".format(runs))\n",
- "print(\"Modules: {}\".format(modules))\n",
- "print(\"Sequences: {}\".format(sequences))\n",
- "print(\"Using DB: {}\".format(db_output))\n",
- "\n",
- "\n",
- "if instrument == \"SPB\":\n",
- " loc = \"SPB_DET_AGIPD1M-1\"\n",
- " dinstance = \"AGIPD1M1\"\n",
- " karabo_id_control = \"SPB_IRU_AGIPD1M1\"\n",
- "else:\n",
- " loc = \"MID_DET_AGIPD1M-1\"\n",
- " dinstance = \"AGIPD1M2\"\n",
- " karabo_id_control = \"MID_EXP_AGIPD1M1\"\n",
- "\n",
- "cal_db_interface = get_random_db_interface(cal_db_interface)\n",
- "\n",
- "# these lines can usually stay as is\n",
- "fbase = \"{}/{{}}/RAW-{{}}-AGIPD{{:02d}}-S{{:05d}}.h5\".format(in_folder)\n",
- "gains = np.arange(3)\n",
- "\n",
- "\n",
- "run, prop, seq = run_prop_seq_from_path(in_folder)\n",
- "logger = InfluxLogger(detector=\"AGIPD\", instrument=instrument, mem_cells=mem_cells,\n",
- " notebook=get_notebook_name(), proposal=prop)\n",
- "\n",
- "# extract the memory cells and acquisition rate from \n",
- "# the file of the first module, first sequence and first run\n",
- "channel = 0\n",
- "fname = fbase.format(runs[0], runs[0].upper(), channel, sequences[0])\n",
- "if acqrate == 0.:\n",
- " acqrate = get_acq_rate(fname, loc, channel)\n",
- " \n",
- "\n",
- "if mem_cells == 0:\n",
- " cells = get_num_cells(fname, loc, channel)\n",
- " mem_cells = cells # avoid setting twice\n",
- " \n",
- "\n",
- "IL_MODE = interlaced\n",
- "max_cells = mem_cells if not interlaced else mem_cells*2\n",
- "memory_cells = mem_cells\n",
- "print(\"Interlaced mode: {}\".format(IL_MODE))\n",
- "\n",
- "cells = np.arange(max_cells)\n",
- "\n",
- "print(f\"Acquisition rate and memory cells set from file {fname} are \" \n",
- " f\"{acqrate} MHz and {max_cells}, respectively\")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "control_fname = f'{in_folder}/{runs[0]}/RAW-{runs[0].upper()}-{karabo_da_control}-S00000.h5'\n",
- "\n",
- "if \"{\" in h5path_ctrl:\n",
- " h5path_ctrl = h5path_ctrl.format(karabo_id_control)\n",
- "\n",
- "if gain_setting == 0.1:\n",
- " if creation_time.replace(tzinfo=None) < dateutil.parser.parse('2020-01-31'):\n",
- " print(\"Set gain-setting to None for runs taken before 2020-01-31\")\n",
- " gain_setting = None\n",
- " else:\n",
- " try:\n",
- " gain_setting = get_gain_setting(control_fname, h5path_ctrl)\n",
- " except Exception as e:\n",
- " print(f'Error while reading gain setting from: \\n{control_fname}')\n",
- " print(e)\n",
- " print(\"Gain setting is not found in the control information\")\n",
- " print(\"Data will not be processed\")\n",
- " sequences = []\n",
- "print(f\"Gain setting: {gain_setting}\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "For the characterization offset maps for each module are needed. The are read in either locally or through the XFEL calibration database."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T07:34:45.876476Z",
- "start_time": "2019-09-13T07:34:25.656979Z"
- }
- },
- "outputs": [],
- "source": [
- "from dateutil import parser\n",
- "offset_g = {}\n",
- "noise_g = {}\n",
- "thresholds_g = {}\n",
- "pc_g = {}\n",
- "if not db_input:\n",
- " print(\"Offset, noise and thresholds have been read in from: {}\".format(offset_store))\n",
- " store_file = h5py.File(offset_store, \"r\")\n",
- " for i in modules:\n",
- " qm = \"Q{}M{}\".format(i//4+1, i%4+1)\n",
- " offset_g[qm] = np.array(store_file[\"{}/Offset/0/data\".format(qm)])\n",
- " noise_g[qm] = np.array(store_file[\"{}/Noise/0/data\".format(qm)])\n",
- " thresholds_g[qm] = np.array(store_file[\"{}/Threshold/0/data\".format(qm)])\n",
- " store_file.close()\n",
- "else:\n",
- " print(\"Offset, noise and thresholds have been read in from calibration database:\")\n",
- " first_ex = True\n",
- " for i in modules:\n",
- " qm = \"Q{}M{}\".format(i//4+1, i%4+1)\n",
- " metadata = ConstantMetaData()\n",
- " offset = Constants.AGIPD.Offset()\n",
- " metadata.calibration_constant = offset\n",
- " det = getattr(Detectors, dinstance)\n",
- "\n",
- " # set the operating condition\n",
- " condition = Conditions.Dark.AGIPD(memory_cells=max_cells, bias_voltage=bias_voltage,\n",
- " acquisition_rate=acqrate, gain_setting=gain_setting)\n",
- " metadata.detector_condition = condition\n",
- "\n",
- " # specify the a version for this constant\n",
- " if creation_time is None:\n",
- " metadata.calibration_constant_version = Versions.Now(device=getattr(det, qm))\n",
- " else:\n",
- " metadata.calibration_constant_version = Versions.Timespan(device=getattr(det, qm), start=creation_time)\n",
- " metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=300000)\n",
- " offset_g[qm] = offset.data\n",
- " if first_ex:\n",
- " print(\"Offset map was injected on: {}\".format(metadata.calibration_constant_version.begin_at))\n",
- " \n",
- " metadata = ConstantMetaData()\n",
- " noise = Constants.AGIPD.Noise()\n",
- " metadata.calibration_constant = noise\n",
- "\n",
- " # set the operating condition\n",
- " condition = Conditions.Dark.AGIPD(memory_cells=max_cells, bias_voltage=bias_voltage,\n",
- " acquisition_rate=acqrate, gain_setting=gain_setting)\n",
- " metadata.detector_condition = condition\n",
- "\n",
- " # specify the a version for this constant\n",
- " if creation_time is None:\n",
- " metadata.calibration_constant_version = Versions.Now(device=getattr(det, qm))\n",
- " else:\n",
- " metadata.calibration_constant_version = Versions.Timespan(device=getattr(det, qm), start=creation_time)\n",
- " metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=300000)\n",
- " noise_g[qm] = noise.data\n",
- " if first_ex:\n",
- " print(\"Noise map was injected on: {}\".format(metadata.calibration_constant_version.begin_at))\n",
- " \n",
- " metadata = ConstantMetaData()\n",
- " thresholds = Constants.AGIPD.ThresholdsDark()\n",
- " metadata.calibration_constant = thresholds\n",
- "\n",
- " # set the operating condition\n",
- " condition = Conditions.Dark.AGIPD(memory_cells=max_cells, bias_voltage=bias_voltage,\n",
- " acquisition_rate=acqrate, gain_setting=gain_setting)\n",
- " metadata.detector_condition = condition\n",
- "\n",
- " # specify the a version for this constant\n",
- " if creation_time is None:\n",
- " metadata.calibration_constant_version = Versions.Now(device=getattr(det, qm))\n",
- " else:\n",
- " metadata.calibration_constant_version = Versions.Timespan(device=getattr(det, qm), start=creation_time)\n",
- " metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=300000)\n",
- " thresholds_g[qm] = thresholds.data\n",
- " if first_ex:\n",
- " print(\"Threshold map was injected on: {}\".format(metadata.calibration_constant_version.begin_at))\n",
- " \n",
- " \n",
- " metadata = ConstantMetaData()\n",
- " pc = Constants.AGIPD.SlopesPC()\n",
- " metadata.calibration_constant = pc\n",
- "\n",
- " # set the operating condition\n",
- " condition = Conditions.Dark.AGIPD(memory_cells=max_cells, bias_voltage=bias_voltage,\n",
- " acquisition_rate=acqrate, gain_setting=gain_setting)\n",
- " metadata.detector_condition = condition\n",
- "\n",
- " # specify the a version for this constant\n",
- " if creation_time is None:\n",
- " metadata.calibration_constant_version = Versions.Now(device=getattr(det, qm))\n",
- " else:\n",
- " metadata.calibration_constant_version = Versions.Timespan(device=getattr(det, qm), start=creation_time)\n",
- " metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=300000)\n",
- " pc_g[qm] = np.nanmedian(pc.data[0,...])/pc.data\n",
- " if first_ex:\n",
- " print(\"PC map was injected on: {}\".format(metadata.calibration_constant_version.begin_at))\n",
- " first_ex = False\n",
- " "
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Initial peak estimates ##\n",
- "\n",
- "First initial peak estimates need to be defined. This is done by inspecting histograms created from (a subset of) the offset-corrected data for peak locations and peak heights."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T07:51:37.815384Z",
- "start_time": "2019-09-13T07:34:45.879121Z"
- }
- },
- "outputs": [],
- "source": [
- "def hist_single_module(fbase, runs, sequences, sensor_size, memory_cells, block_size,\n",
- " il_mode, limit_trains, profile, rawversion, instrument, inp):\n",
- " \"\"\" This function calculates a per-pixel histogram for a single module\n",
- " \n",
- " Runs and sequences give the data to calculate histogram from\n",
- " \"\"\"\n",
- " channel, offset, thresholds, pc, noise = inp\n",
- " \n",
- " import XFELDetAna.xfelpycaltools as xcal\n",
- " import numpy as np\n",
- " import h5py\n",
- " from XFELDetAna.util import env\n",
- " \n",
- " \n",
- " env.iprofile = profile\n",
- " \n",
- " def baseline_correct_via_noise(d, noise, g, baseline_corr_noise_threshold=1000):\n",
- " \"\"\" Correct baseline shifts by evaluating position of the noise peak\n",
- "\n",
- " :param: d: the data to correct, should be a single image\n",
- " :param noise: the mean noise for the cell id of `d`\n",
- " :param g: gain array matching `d` array\n",
- "\n",
- " Correction is done by identifying the left-most (significant) peak\n",
- " in the histogram of `d` and shifting it to be centered on 0.\n",
- " This is done via a continous wavelet transform (CWT), using a Ricker\n",
- " wavelet.\n",
- "\n",
- " Only high gain data is evaulated, and data larger than 50 ADU,\n",
- " equivalent of roughly a 9 keV photon is ignored.\n",
- "\n",
- " Two passes are executed: the first shift is accurate to 10 ADU and\n",
- " will catch baseline shifts smaller then -2000 ADU. CWT is evaluated\n",
- " for peaks of widths one, three and five time the noise.\n",
- " The baseline is then shifted by the position of the summed CWT maximum.\n",
- "\n",
- " In a second pass histogram is evaluated within a range of\n",
- " +- 5*noise of the initial shift value, for peak of width noise.\n",
- "\n",
- " Baseline shifts larger than the maximum threshold or positive shifts\n",
- " larger 10 are discarded, and the original data is returned, otherwise\n",
- " the shift corrected data is returned.\n",
- "\n",
- " \"\"\"\n",
- " import copy\n",
- " from scipy.signal import cwt, ricker\n",
- " # we work on a copy to be able to filter values by setting them to\n",
- " # np.nan\n",
- " dd = copy.copy(d)\n",
- " dd[g != 0] = np.nan # only high gain data\n",
- " dd[dd > 50] = np.nan # only noise data\n",
- " h, e = np.histogram(dd, bins=210, range=(-2000, 100))\n",
- " c = (e[1:] + e[:-1]) / 2\n",
- "\n",
- " try:\n",
- " cwtmatr = cwt(h, ricker, [noise, 3. * noise, 5. * noise])\n",
- " except:\n",
- " return d\n",
- " cwtall = np.sum(cwtmatr, axis=0)\n",
- " marg = np.argmax(cwtall)\n",
- " pc = c[marg]\n",
- "\n",
- " high_res_range = (dd > pc - 5 * noise) & (dd < pc + 5 * noise)\n",
- " dd[~high_res_range] = np.nan\n",
- " h, e = np.histogram(dd, bins=200,\n",
- " range=(pc - 5 * noise, pc + 5 * noise))\n",
- " c = (e[1:] + e[:-1]) / 2\n",
- " try:\n",
- " cwtmatr = cwt(h, ricker, [noise, ])\n",
- " except:\n",
- " return d\n",
- " marg = np.argmax(cwtmatr)\n",
- " pc = c[marg]\n",
- " # too large shift to be easily decernable via noise\n",
- " if pc > 10 or pc < -baseline_corr_noise_threshold:\n",
- " return d\n",
- " return d - pc\n",
- "\n",
- " \n",
- " # function needs to be inline for parallell processing\n",
- " def read_fun(filename, channel):\n",
- " \"\"\" A reader function used by pyDetLib\n",
- " \"\"\"\n",
- " infile = h5py.File(filename, \"r\", driver=\"core\")\n",
- " if rawversion == 2:\n",
- " count = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/count\".format(instrument, channel)])\n",
- " first = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/first\".format(instrument, channel)])\n",
- " last_index = int(first[count != 0][-1]+count[count != 0][-1])\n",
- " first_index = int(first[count != 0][0])\n",
- " else:\n",
- " status = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/status\".format(instrument, channel)])\n",
- " if np.count_nonzero(status != 0) == 0:\n",
- " return\n",
- " last = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/last\".format(instrument, channel)])\n",
- " last_index = int(last[status != 0][-1])\n",
- " first_index = int(last[status != 0][0])\n",
- " if limit_trains is not None:\n",
- " last_index = min(limit_trains*memory_cells+first_index, last_index)\n",
- " \n",
- " im = np.array(infile[\"/INSTRUMENT/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/data\".format(instrument, channel)][first_index:last_index,...]) \n",
- " carr = infile[\"/INSTRUMENT/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/cellId\".format(instrument, channel)][first_index:last_index]\n",
- " cells = np.squeeze(np.array(carr))\n",
- " infile.close()\n",
- " \n",
- " if il_mode:\n",
- " ga = im[1::2, 0, ...]\n",
- " im = im[0::2, 0, ...].astype(np.float32)\n",
- " else:\n",
- " ga = im[:, 1, ...]\n",
- " im = im[:, 0, ...].astype(np.float32)\n",
- "\n",
- " im = np.rollaxis(im, 2)\n",
- " im = np.rollaxis(im, 2, 1)\n",
- "\n",
- " ga = np.rollaxis(ga, 2)\n",
- " ga = np.rollaxis(ga, 2, 1)\n",
- " return im, ga, cells\n",
- " \n",
- " offset_cor = xcal.OffsetCorrection(sensor_size,\n",
- " offset,\n",
- " nCells=memory_cells,\n",
- " blockSize=block_size,\n",
- " gains=[0,1,2])\n",
- " offset_cor.mapper = offset_cor._view.map_sync\n",
- " offset_cor.debug() # force non-parallel processing since outer function will run concurrently\n",
- " hist_calc = xcal.HistogramCalculator(sensor_size,\n",
- " bins=4000,\n",
- " range=(-4000, 8000),\n",
- " blockSize=block_size)\n",
- " hist_calc.mapper = hist_calc._view.map_sync\n",
- " hist_calc.debug() # force non-parallel processing since outer function will run concurrently\n",
- " for run in runs:\n",
- " for seq in sequences:\n",
- " fname = fbase.format(run, run.upper(), channel, seq)\n",
- " d, ga, c = read_fun(fname, channel)\n",
- " vidx = (c < offset.shape[2]) & (c < thresholds.shape[2])\n",
- " c = c[vidx]\n",
- " d = d[:,:,vidx]\n",
- " ga = ga[:,:,vidx]\n",
- " # we need to do proper gain thresholding\n",
- " g = np.zeros(ga.shape, np.uint8)\n",
- " g[...] = 2\n",
- " \n",
- " g[ga < thresholds[...,c,1]] = 1\n",
- " g[ga < thresholds[...,c,0]] = 0\n",
- " d = offset_cor.correct(d, cellTable=c, gainMap=g)\n",
- " for i in range(d.shape[2]):\n",
- " mn_noise = np.nanmean(noise[...,c[i]])\n",
- " d[...,i] = baseline_correct_via_noise(d[...,i],\n",
- " mn_noise,\n",
- " g[..., i])\n",
- " \n",
- " d *= np.moveaxis(pc[c,...], 0, 2)\n",
- " \n",
- " hist_calc.fill(d)\n",
- " h, e, c, _ = hist_calc.get()\n",
- " return h, e, c\n",
- "\n",
- "inp = []\n",
- "for i in modules:\n",
- " qm = \"Q{}M{}\".format(i//4+1, i%4+1)\n",
- " inp.append((i, offset_g[qm], thresholds_g[qm], pc_g[qm][0,...], noise_g[qm][...,0]))\n",
- " \n",
- "p = partial(hist_single_module, fbase, runs, sequences,\n",
- " sensor_size, memory_cells, block_size, IL_MODE, limit_trains, profile, rawversion, instrument)\n",
- "res_uncorr = list(map(p, inp))\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T07:51:38.024035Z",
- "start_time": "2019-09-13T07:51:37.818276Z"
- },
- "scrolled": false
- },
- "outputs": [],
- "source": [
- "d = []\n",
- "qms = []\n",
- "for i, r in enumerate(res_uncorr):\n",
- " ii = list(modules)[i]\n",
- " qm = \"Q{}M{}\".format(ii//4+1, ii%4+1)\n",
- " qms.append(qm)\n",
- " h, e, c = r \n",
- " d.append({\n",
- " 'x': c, \n",
- " 'y': h,\n",
- " 'drawstyle': 'steps-mid'\n",
- " })\n",
- " \n",
- "fig = xana.simplePlot(d, y_log=False,\n",
- " figsize=\"2col\",\n",
- " aspect=2,\n",
- " x_range=(-50, 500),\n",
- " x_label=\"Intensity (ADU)\",\n",
- " y_label=\"Counts\")\n",
- "\n",
- "fig.savefig(\"{}/FF_module_{}_peak_pos.png\".format(out_folder, \",\".join(qms)))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T07:51:38.029704Z",
- "start_time": "2019-09-13T07:51:38.026111Z"
- }
- },
- "outputs": [],
- "source": [
- "# these should be quite stable\n",
- "\n",
- "peak_estimates = [0, 55, 105, 155]\n",
- "peak_heights = [5e7, 5e6, 1e6, 5e5]\n",
- "peak_sigma = [5., 5.,5., 5.,]\n",
- "\n",
- "print(\"Using the following peak estimates: {}\".format(peak_estimates))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "collapsed": true
- },
- "source": [
- "## Calculate relative gain per module ##\n",
- "\n",
- "Using the so obtained estimates, we calculate the relative gain per module. For this we use the weighted average method implemented in pyDetLib."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T08:16:50.652492Z",
- "start_time": "2019-09-13T07:51:38.031787Z"
- }
- },
- "outputs": [],
- "source": [
- "block_size = [64, 64]\n",
- "def relgain_single_module(fbase, runs, sequences, peak_estimates,\n",
- " peak_heights, peak_sigma, memory_cells, sensor_size,\n",
- " block_size, il_mode, profile, limit_trains_eval, rawversion, instrument, inp):\n",
- " \"\"\" A function for calculated the relative gain of a single AGIPD module\n",
- " \"\"\"\n",
- " \n",
- " # import needed inline for parallel processing\n",
- " import XFELDetAna.xfelpycaltools as xcal\n",
- " import numpy as np\n",
- " import h5py\n",
- " from XFELDetAna.util import env\n",
- " env.iprofile = profile\n",
- " \n",
- " channel, offset, thresholds, noise, pc = inp\n",
- " \n",
- " def baseline_correct_via_noise(d, noise, g, baseline_corr_noise_threshold=1000):\n",
- " \"\"\" Correct baseline shifts by evaluating position of the noise peak\n",
- "\n",
- " :param: d: the data to correct, should be a single image\n",
- " :param noise: the mean noise for the cell id of `d`\n",
- " :param g: gain array matching `d` array\n",
- "\n",
- " Correction is done by identifying the left-most (significant) peak\n",
- " in the histogram of `d` and shifting it to be centered on 0.\n",
- " This is done via a continous wavelet transform (CWT), using a Ricker\n",
- " wavelet.\n",
- "\n",
- " Only high gain data is evaulated, and data larger than 50 ADU,\n",
- " equivalent of roughly a 9 keV photon is ignored.\n",
- "\n",
- " Two passes are executed: the first shift is accurate to 10 ADU and\n",
- " will catch baseline shifts smaller then -2000 ADU. CWT is evaluated\n",
- " for peaks of widths one, three and five time the noise.\n",
- " The baseline is then shifted by the position of the summed CWT maximum.\n",
- "\n",
- " In a second pass histogram is evaluated within a range of\n",
- " +- 5*noise of the initial shift value, for peak of width noise.\n",
- "\n",
- " Baseline shifts larger than the maximum threshold or positive shifts\n",
- " larger 10 are discarded, and the original data is returned, otherwise\n",
- " the shift corrected data is returned.\n",
- "\n",
- " \"\"\"\n",
- " import copy\n",
- " from scipy.signal import cwt, ricker\n",
- " # we work on a copy to be able to filter values by setting them to\n",
- " # np.nan\n",
- " dd = copy.copy(d)\n",
- " dd[g != 0] = np.nan # only high gain data\n",
- " dd[dd > 50] = np.nan # only noise data\n",
- " h, e = np.histogram(dd, bins=210, range=(-2000, 100))\n",
- " c = (e[1:] + e[:-1]) / 2\n",
- "\n",
- " try:\n",
- " cwtmatr = cwt(h, ricker, [noise, 3. * noise, 5. * noise])\n",
- " except:\n",
- " return d\n",
- " cwtall = np.sum(cwtmatr, axis=0)\n",
- " marg = np.argmax(cwtall)\n",
- " pc = c[marg]\n",
- "\n",
- " high_res_range = (dd > pc - 5 * noise) & (dd < pc + 5 * noise)\n",
- " dd[~high_res_range] = np.nan\n",
- " h, e = np.histogram(dd, bins=200,\n",
- " range=(pc - 5 * noise, pc + 5 * noise))\n",
- " c = (e[1:] + e[:-1]) / 2\n",
- " try:\n",
- " cwtmatr = cwt(h, ricker, [noise, ])\n",
- " except:\n",
- " return d\n",
- " marg = np.argmax(cwtmatr)\n",
- " pc = c[marg]\n",
- " # too large shift to be easily decernable via noise\n",
- " if pc > 10 or pc < -baseline_corr_noise_threshold:\n",
- " return d\n",
- " return d - pc\n",
- "\n",
- " \n",
- " # function needs to be inline for parallell processing\n",
- " def read_fun(filename, channel):\n",
- " \"\"\" A reader function used by pyDetLib\n",
- " \"\"\"\n",
- " infile = h5py.File(filename, \"r\", driver=\"core\")\n",
- " if rawversion == 2:\n",
- " count = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/count\".format(instrument, channel)])\n",
- " first = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/first\".format(instrument, channel)])\n",
- " last_index = int(first[count != 0][-1]+count[count != 0][-1])\n",
- " first_index = int(first[count != 0][0])\n",
- " else:\n",
- " status = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/status\".format(instrument, channel)])\n",
- " if np.count_nonzero(status != 0) == 0:\n",
- " return\n",
- " last = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/last\".format(instrument, channel)])\n",
- " last_index = int(last[status != 0][-1])\n",
- " first_index = int(last[status != 0][0])\n",
- " if limit_trains is not None:\n",
- " last_index = min(limit_trains*memory_cells+first_index, last_index)\n",
- " im = np.array(infile[\"/INSTRUMENT/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/data\".format(instrument, channel)][first_index:last_index,...]) \n",
- " carr = infile[\"/INSTRUMENT/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/cellId\".format(instrument, channel)][first_index:last_index]\n",
- " cells = np.squeeze(np.array(carr))\n",
- " infile.close()\n",
- "\n",
- " if il_mode:\n",
- " ga = im[1::2, 0, ...]\n",
- " im = im[0::2, 0, ...].astype(np.float32)\n",
- " else:\n",
- " ga = im[:, 1, ...]\n",
- " im = im[:, 0, ...].astype(np.float32)\n",
- " \n",
- " im = np.rollaxis(im, 2)\n",
- " im = np.rollaxis(im, 2, 1)\n",
- "\n",
- " ga = np.rollaxis(ga, 2)\n",
- " ga = np.rollaxis(ga, 2, 1)\n",
- " return im, ga, cells\n",
- " \n",
- " offset_cor = xcal.OffsetCorrection(sensor_size, offset, nCells=memory_cells,\n",
- " blockSize=block_size, gains=[0,1,2])\n",
- " offset_cor.mapper = offset_cor._view.map_sync\n",
- "\n",
- " rel_gain = xcal.GainMapCalculator(sensor_size,\n",
- " peak_estimates,\n",
- " peak_sigma,\n",
- " nCells=memory_cells,\n",
- " peakHeights = peak_heights,\n",
- " noiseMap=noise,\n",
- " deviationType=\"relative\")\n",
- " rel_gain.mapper = rel_gain._view.map_sync\n",
- " for run in runs:\n",
- " for seq in sequences:\n",
- " fname = fbase.format(run, run.upper(), channel, seq)\n",
- " d, ga, c = read_fun(fname, channel)\n",
- " # we need to do proper gain thresholding\n",
- " vidx = (c < offset.shape[2]) & (c < thresholds.shape[2])\n",
- " c = c[vidx]\n",
- " d = d[:,:,vidx]\n",
- " ga = ga[:,:,vidx]\n",
- " # we need to do proper gain thresholding\n",
- " g = np.zeros(ga.shape, np.uint8)\n",
- " g[...] = 2\n",
- " \n",
- " g[ga < thresholds[...,c,1]] = 1\n",
- " g[ga < thresholds[...,c,0]] = 0\n",
- " d = offset_cor.correct(d, cellTable=c, gainMap=g)\n",
- " \n",
- " for i in range(d.shape[2]):\n",
- " mn_noise = np.nanmean(noise[...,c[i]])\n",
- " d[...,i] = baseline_correct_via_noise(d[...,i],\n",
- " mn_noise,\n",
- " g[..., i])\n",
- " \n",
- " d *= np.moveaxis(pc[c,...], 0, 2)\n",
- " rel_gain.fill(d, cellTable=c)\n",
- " \n",
- " gain_map = rel_gain.get()\n",
- " gain_map_func = rel_gain.getUsingFunc(inverse=False)\n",
- " \n",
- " pks, stds = rel_gain.getRawPeaks()\n",
- " return gain_map, pks, stds, gain_map_func\n",
- "\n",
- "inp = []\n",
- "for i in modules:\n",
- " qm = \"Q{}M{}\".format(i//4+1, i%4+1)\n",
- " inp.append((i, offset_g[qm], thresholds_g[qm], noise_g[qm][...,0], pc_g[qm][0,...]))\n",
- "start = datetime.now()\n",
- "p = partial(relgain_single_module, fbase, runs, sequences,\n",
- " peak_estimates, peak_heights, peak_sigma, memory_cells,\n",
- " sensor_size, block_size, IL_MODE, profile, limit_trains_eval, rawversion, instrument)\n",
- "res_gain = list(map(p, inp)) # don't run concurently as inner function are parelllized\n",
- "duration = (datetime.now()-start).total_seconds()\n",
- "logger.runtime_summary_entry(success=True, runtime=duration)\n",
- "logger.send()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Finally, we inspect the results, by plotting the number of entries, peak locations and resulting gain maps for each peak. In the course of doing so, we identify bad pixels by either having 0 entries for a peak, or having `nan` values for the peak location."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T08:16:54.261782Z",
- "start_time": "2019-09-13T08:16:50.657594Z"
- }
- },
- "outputs": [],
- "source": [
- "\n",
- "gain_m = {}\n",
- "flatsff = {}\n",
- "gainoff_g = {}\n",
- "entries_g = {}\n",
- "peaks_g = {}\n",
- "mask_g = {}\n",
- "cell_to_preview = 4\n",
- "masks_eval_peaks = [1, 2]\n",
- "global_eval_peaks = [1]\n",
- "global_meds = {}\n",
- "\n",
- "for i, r in enumerate(res_gain):\n",
- " ii = list(modules)[i]\n",
- " qm = \"Q{}M{}\".format(ii//4+1, ii%4+1)\n",
- " gain, pks, std, gfunc = r\n",
- " gains, errors, chisq, valid, max_dev, stats = gfunc\n",
- " _, entries, stds, sow = gain\n",
- " gain_db = np.zeros((gains.shape[0], gains.shape[1], memory_cells))\n",
- " gain_mdb = np.zeros((gains.shape[0], gains.shape[1], memory_cells))\n",
- " entries_db = np.zeros((gains.shape[0], gains.shape[1], memory_cells, 5)) \n",
- " gainoff_db = np.zeros((gains.shape[0], gains.shape[1], memory_cells))\n",
- " mask_db = np.zeros((gains.shape[0], gains.shape[1], memory_cells), np.uint32)\n",
- " \n",
- " gainoff_g[qm] = gainoff_db\n",
- " gain_m[qm] = gain_mdb\n",
- " entries_g[qm] = entries_db\n",
- " peaks_g[qm] = pks\n",
- " \n",
- " # create a mask for unregular pixels\n",
- " # first bit set if first peak has nan entries\n",
- " for pk in masks_eval_peaks:\n",
- " mask_db[~(np.isfinite(gain_mdb) & np.isfinite(gainoff_db))] |= BadPixels.FF_GAIN_EVAL_ERROR.value\n",
- " mask_db[(np.abs(1-gain_mdb/np.nanmedian(gain_mdb) > deviation_threshold) )] |= BadPixels.FF_GAIN_DEVIATION.value\n",
- " # second bit set if entries are 0 for first peak\n",
- " mask_db[entries[...,1] == 0] |= BadPixels.FF_NO_ENTRIES.value\n",
- " zero_oc = np.count_nonzero((entries > 0).astype(np.int), axis=3)\n",
- " mask_db[zero_oc <= len(peak_estimates)-2] |= BadPixels.FF_NO_ENTRIES.value \n",
- " # third bit set if entries of a given adc show significant noise\n",
- " stds = []\n",
- " for ii in range(8):\n",
- " for jj in range(8):\n",
- " stds.append(np.std(entries_db[ii*16:(ii+1)*16,jj*64+2:(jj+1)*64-2,:,1], axis=(0,1)))\n",
- " avg_stds = np.median(np.array(stds), axis=0)\n",
- " \n",
- " for ii in range(8):\n",
- " for jj in range(8):\n",
- " std = np.std(entries_db[ii*16:(ii+1)*16,jj*64+2:(jj+1)*64-2,:,1], axis=(0,1))\n",
- " if np.any(std > 5*avg_stds):\n",
- " mask_db[ii*16:(ii+1)*16,jj*64:(jj+1)*64,std > 5*avg_stds] |= BadPixels.FF_GAIN_DEVIATION\n",
- " \n",
- " \n",
- " mask_g[qm] = mask_db\n",
- " \n",
- " flat = np.zeros((gains.shape[0], gains.shape[1], memory_cells, 3))\n",
- " for j in range(2,len(peak_estimates)):\n",
- " flat[...,j-2] = np.mean(entries[...,j]/np.mean(entries[...,j]))\n",
- " flat = np.mean(flat, axis=3)\n",
- " #flat_db = np.zeros((gains.shape[0], gains.shape[1], memory_cells))\n",
- " #for j in range(2):\n",
- " # flat_db[...,j::2] = flat\n",
- " flatsff[qm] = flat\n",
- " \n",
- " global_meds[qm] = np.nanmedian(pks[...,global_eval_peaks][np.max(mask_db, axis=2) != 0])"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Evaluated peak locations ##\n",
- "\n",
- "The following plot shows the evaluated peak locations for each peak. Peak ids increase downward:\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T08:16:55.562017Z",
- "start_time": "2019-09-13T08:16:54.264843Z"
- },
- "scrolled": false
- },
- "outputs": [],
- "source": [
- "from mpl_toolkits.axes_grid1 import AxesGrid\n",
- "cell_to_preview = 4\n",
- "for qm, data in peaks_g.items():\n",
- " print(\"The module shown is: {}\".format(qm))\n",
- " print(\"The cell shown is: {}\".format(cell_to_preview))\n",
- " print(\"Evaluated peaks: {}\".format(data.shape[-1]))\n",
- " fig = plt.figure(figsize=(15,15))\n",
- " grid = AxesGrid(fig, 111, \n",
- " nrows_ncols=(data.shape[-1], 1),\n",
- " axes_pad=0.1,\n",
- " share_all=True,\n",
- " label_mode=\"L\",\n",
- " cbar_location=\"right\",\n",
- " cbar_mode=\"each\",\n",
- " cbar_size=\"7%\",\n",
- " cbar_pad=\"2%\")\n",
- " \n",
- " for j in range(data.shape[-1]):\n",
- " d = data[...,cell_to_preview,j]\n",
- " d[~np.isfinite(d)] = 0\n",
- " \n",
- " im = grid[j].imshow(d, interpolation=\"nearest\", vmin=0.9*np.nanmedian(d), vmax=max(1.1*np.nanmedian(d), 50))\n",
- " cb = grid.cbar_axes[j].colorbar(im)\n",
- " cb.set_label_text(\"Peak location (ADU)\")"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Gain Slopes And Offsets ##\n",
- "\n",
- "The gain slopes and offsets are deduced by fitting a linear function $y = mx + b$ to the evaluated peaks. Gains are normalized to all pixels and all memory cells of a module by using the average peak locations a $x$ values. Thus slopes centered around 1 are to be expected.\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T08:16:56.879582Z",
- "start_time": "2019-09-13T08:16:55.564572Z"
- },
- "scrolled": false
- },
- "outputs": [],
- "source": [
- "for i, r in enumerate(res_gain):\n",
- " ii = list(modules)[i]\n",
- " qm = \"Q{}M{}\".format(ii//4+1, ii%4+1)\n",
- " gain, pks, std, gfunc = r\n",
- " gains, errors, chisq, valid, max_dev, stats = gfunc \n",
- " _, entries, stds, sow = gain\n",
- " \n",
- " fig = plt.figure(figsize=(15,8))\n",
- " ax = fig.add_subplot(211)\n",
- " im = ax.imshow(gains[...,0], interpolation=\"nearest\", vmin=0.85, vmax=1.15)\n",
- " cb = fig.colorbar(im)\n",
- " cb.set_label(\"Gain slope m\")\n",
- " fig.savefig(\"{}/gain_m_mod{}.png\".format(out_folder, qm))\n",
- " \n",
- " ax = fig.add_subplot(212)\n",
- " ax.hist(gains[...,0].flatten(), range=(0.5, 1.5), bins=100)\n",
- " ax.set_ylabel(\"Occurences\")\n",
- " ax.set_xlabel(\"Gain slope m\")\n",
- " ax.semilogy()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The gain offsets b are expected to be centered around 0 and minimal, as data is already offset substracted."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T08:16:58.113910Z",
- "start_time": "2019-09-13T08:16:56.881383Z"
- }
- },
- "outputs": [],
- "source": [
- "for i, r in enumerate(res_gain):\n",
- " ii = list(modules)[i]\n",
- " qm = \"Q{}M{}\".format(ii//4+1, ii%4+1)\n",
- " gain, pks, std, gfunc = r\n",
- " gains, errors, chisq, valid, max_dev, stats = gfunc\n",
- " _, entries, stds, sow = gain\n",
- " \n",
- " fig = plt.figure(figsize=(15,8))\n",
- " ax = fig.add_subplot(211)\n",
- " \n",
- " im = ax.imshow(gains[...,1], interpolation=\"nearest\", vmin=-25, vmax=25)\n",
- " cb = fig.colorbar(im)\n",
- " cb.set_label(\"Gain offset b\")\n",
- " fig.savefig(\"{}/gain_b_mod{}.png\".format(out_folder, qm))\n",
- " \n",
- " ax = fig.add_subplot(212)\n",
- " ax.hist(gains[...,1].flatten(), range=(-25, 25), bins=100)\n",
- " ax.set_ylabel(\"Occurences\")\n",
- " ax.set_xlabel(\"Gain offset b\")\n",
- " ax.semilogy()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Bad Pixels ##\n",
- "\n",
- "Bad pixels are defined as any of the following criteria: \n",
- "\n",
- "* the gain evaluation did not converge, or location of noise peak deviates by more than than the deviation threshold from the median location.\n",
- "* the number of entries for the noise peak in 0.\n",
- "* the standard deviation of the number of entries is larger than 5 times the standard deviation for the ASIC the pixel is on."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T08:18:08.277881Z",
- "start_time": "2019-09-13T08:18:08.272390Z"
- }
- },
- "outputs": [],
- "source": [
- "print(\"The deviation threshold is: {}\".format(deviation_threshold))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T08:18:10.840520Z",
- "start_time": "2019-09-13T08:18:09.769451Z"
- }
- },
- "outputs": [],
- "source": [
- "for i, r in enumerate(res_gain):\n",
- " ii = list(modules)[i]\n",
- " qm = \"Q{}M{}\".format(ii//4+1, ii%4+1)\n",
- " mask_db = mask_g[qm]\n",
- " \n",
- " fig = plt.figure(figsize=(15,8))\n",
- " ax = fig.add_subplot(211)\n",
- " im = ax.imshow(np.log2(mask_db[...,cell_to_preview]), interpolation=\"nearest\", vmin=0, vmax=32)\n",
- " cb = fig.colorbar(im)\n",
- " cb.set_label(\"Mask value\")\n",
- " fig.savefig(\"{}/mask_mod{}.png\".format(out_folder, qm))\n",
- " \n",
- " ax = fig.add_subplot(212)\n",
- " ax.hist(np.log2(mask_db.flatten()), range=(0, 32), bins=32, normed=True)\n",
- " ax.set_ylabel(\"Occurences\")\n",
- " ax.set_xlabel(\"Mask value\")\n",
- " ax.semilogy()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T08:18:14.349244Z",
- "start_time": "2019-09-13T08:18:11.849053Z"
- }
- },
- "outputs": [],
- "source": [
- "cols = {BadPixels.FF_GAIN_EVAL_ERROR.value: (BadPixels.FF_GAIN_EVAL_ERROR.name, '#FF000080'),\n",
- " BadPixels.FF_GAIN_DEVIATION.value: (BadPixels.FF_GAIN_DEVIATION.name, '#0000FF80'),\n",
- " BadPixels.FF_NO_ENTRIES.value: (BadPixels.FF_NO_ENTRIES.name, '#00FF0080'),\n",
- " BadPixels.FF_GAIN_EVAL_ERROR.value |\n",
- " BadPixels.FF_GAIN_DEVIATION.value: ('EVAL+DEV', '#DD00DD80'),\n",
- " BadPixels.FF_GAIN_EVAL_ERROR.value |\n",
- " BadPixels.FF_NO_ENTRIES.value: ('EVAL+NO_ENTRY', '#00DDDD80'),\n",
- " BadPixels.FF_GAIN_DEVIATION.value |\n",
- " BadPixels.FF_NO_ENTRIES.value: ('DEV+NO_ENTRY', '#DDDD0080'),\n",
- " }\n",
- "\n",
- "rebin = 32 if not high_res_badpix_3d else 2\n",
- "\n",
- "gain = 0\n",
- "for mod, data in mask_g.items():\n",
- " plot_badpix_3d(data, cols, title=mod, rebin_fac=rebin)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-07-23T18:02:58.599579Z",
- "start_time": "2019-07-23T18:02:58.316311Z"
- }
- },
- "outputs": [],
- "source": [
- "if local_output:\n",
- " ofile = \"{}/agipd_gain_store_{}_modules_{}.h5\".format(out_folder,\n",
- " \"_\".join([str(r) for r in runs]),\n",
- " \"_\".join([str(m) for m in modules]))\n",
- " store_file = h5py.File(ofile, \"w\")\n",
- " for i, r in enumerate(res_gain):\n",
- " ii = list(modules)[i]\n",
- " qm = \"Q{}M{}\".format(ii//4+1, ii%4+1)\n",
- " gain, pks, std, gfunc = r\n",
- " gains, errors, chisq, valid, max_dev, stats = gfunc\n",
- " gainmap, entires, stds, sow = gain\n",
- " store_file[\"/{}/Gain/0/data\".format(qm)] = gains[...,0]\n",
- " store_file[\"/{}/GainOffset/0/data\".format(qm)] = gains[...,1]\n",
- " store_file[\"/{}/Flat/0/data\".format(qm)] = flatsff[qm]\n",
- " store_file[\"/{}/Entries/0/data\".format(qm)] = entires\n",
- " store_file[\"/{}/BadPixels/0/data\".format(qm)] = mask_g[qm] \n",
- " store_file.close()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]\n",
- "file_loc = proposal + ' ' + ' '.join(list(map(str,runs)))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-12T14:49:28.355922Z",
- "start_time": "2019-09-12T14:49:28.340063Z"
- }
- },
- "outputs": [],
- "source": [
- "if db_output:\n",
- " for i, r in enumerate(res_gain):\n",
- " ii = list(modules)[i]\n",
- " qm = \"Q{}M{}\".format(ii//4+1, ii%4+1)\n",
- " \n",
- " gain, pks, std, gfunc = r\n",
- " gains, errors, chisq, valid, max_dev, stats = gfunc\n",
- " gainmap, entires, stds, sow = gain\n",
- " \n",
- " det = getattr(Detectors, dinstance)\n",
- " device = getattr(det, qm)\n",
- " # gains related\n",
- " metadata = ConstantMetaData()\n",
- " cgain = Constants.AGIPD.SlopesFF()\n",
- " cgain.data = gains\n",
- " metadata.calibration_constant = cgain\n",
- "\n",
- " # set the operating condition\n",
- " condition = Conditions.Illuminated.AGIPD(memory_cells, bias_voltage, 9.2,\n",
- " pixels_x=512, pixels_y=128, beam_energy=None,\n",
- " acquisition_rate=acqrate, gain_setting=gain_setting)\n",
- " \n",
- " metadata.detector_condition = condition\n",
- "\n",
- " # specify the a version for this constant\n",
- " if creation_time is None:\n",
- " metadata.calibration_constant_version = Versions.Now(device=getattr(det, qm))\n",
- " else:\n",
- " metadata.calibration_constant_version = Versions.Timespan(device=getattr(det, qm), start=creation_time)\n",
- "\n",
- " metadata.calibration_constant_version.raw_data_location = file_loc\n",
- " metadata.send(cal_db_interface, timeout=300000)\n",
- " \n",
- " \n",
- " # bad pixels \n",
- " metadata = ConstantMetaData()\n",
- " bp = Constants.AGIPD.BadPixelsFF()\n",
- " bp.data = mask_g[qm] \n",
- " metadata.calibration_constant = bp\n",
- "\n",
- " # set the operating condition\n",
- " condition = Conditions.Illuminated.AGIPD(memory_cells, bias_voltage, 9.2,\n",
- " pixels_x=512, pixels_y=128, beam_energy=None,\n",
- " acquisition_rate=acqrate, gain_setting=gain_setting)\n",
- " \n",
- " metadata.detector_condition = condition\n",
- "\n",
- " # specify the a version for this constant\n",
- " if creation_time is None:\n",
- " metadata.calibration_constant_version = Versions.Now(device=getattr(det, qm))\n",
- " else:\n",
- " metadata.calibration_constant_version = Versions.Timespan(device=getattr(det, qm), start=creation_time)\n",
- " metadata.calibration_constant_version.raw_data_location = file_loc\n",
- " metadata.send(cal_db_interface, timeout=300000)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Sanity check ##\n",
- "\n",
- "Finally, we perform a correction of the data used to derive the gain constants with said constants. We expected an improvement in peak FWHM, if constants have been deduced correctly."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T09:23:54.999086Z",
- "start_time": "2019-09-13T09:16:57.293565Z"
- },
- "scrolled": false
- },
- "outputs": [],
- "source": [
- "def hist_single_module(fbase, runs, sequences, il_mode, profile, limit_trains, memory_cells, rawversion, instrument, inp):\n",
- " channel, offset, thresholds, relgain, noise, pc = inp\n",
- " gain, pks, std, gfunc = relgain\n",
- " gains, errors, chisq, valid, max_dev, stats = gfunc\n",
- " \n",
- " import XFELDetAna.xfelpycaltools as xcal\n",
- " import numpy as np\n",
- " import h5py\n",
- " import copy\n",
- " from XFELDetAna.util import env\n",
- " env.iprofile = profile\n",
- " sensor_size = [128, 512]\n",
- " \n",
- " block_size = sensor_size\n",
- " \n",
- " \n",
- " def baseline_correct_via_noise(d, noise, g, baseline_corr_noise_threshold=1000):\n",
- " \"\"\" Correct baseline shifts by evaluating position of the noise peak\n",
- "\n",
- " :param: d: the data to correct, should be a single image\n",
- " :param noise: the mean noise for the cell id of `d`\n",
- " :param g: gain array matching `d` array\n",
- "\n",
- " Correction is done by identifying the left-most (significant) peak\n",
- " in the histogram of `d` and shifting it to be centered on 0.\n",
- " This is done via a continous wavelet transform (CWT), using a Ricker\n",
- " wavelet.\n",
- "\n",
- " Only high gain data is evaulated, and data larger than 50 ADU,\n",
- " equivalent of roughly a 9 keV photon is ignored.\n",
- "\n",
- " Two passes are executed: the first shift is accurate to 10 ADU and\n",
- " will catch baseline shifts smaller then -2000 ADU. CWT is evaluated\n",
- " for peaks of widths one, three and five time the noise.\n",
- " The baseline is then shifted by the position of the summed CWT maximum.\n",
- "\n",
- " In a second pass histogram is evaluated within a range of\n",
- " +- 5*noise of the initial shift value, for peak of width noise.\n",
- "\n",
- " Baseline shifts larger than the maximum threshold or positive shifts\n",
- " larger 10 are discarded, and the original data is returned, otherwise\n",
- " the shift corrected data is returned.\n",
- "\n",
- " \"\"\"\n",
- " import copy\n",
- " from scipy.signal import cwt, ricker\n",
- " # we work on a copy to be able to filter values by setting them to\n",
- " # np.nan\n",
- " dd = copy.copy(d)\n",
- " dd[g != 0] = np.nan # only high gain data\n",
- " dd[dd > 50] = np.nan # only noise data\n",
- " h, e = np.histogram(dd, bins=210, range=(-2000, 100))\n",
- " c = (e[1:] + e[:-1]) / 2\n",
- "\n",
- " try:\n",
- " cwtmatr = cwt(h, ricker, [noise, 3. * noise, 5. * noise])\n",
- " except:\n",
- " return d\n",
- " cwtall = np.sum(cwtmatr, axis=0)\n",
- " marg = np.argmax(cwtall)\n",
- " pc = c[marg]\n",
- "\n",
- " high_res_range = (dd > pc - 5 * noise) & (dd < pc + 5 * noise)\n",
- " dd[~high_res_range] = np.nan\n",
- " h, e = np.histogram(dd, bins=200,\n",
- " range=(pc - 5 * noise, pc + 5 * noise))\n",
- " c = (e[1:] + e[:-1]) / 2\n",
- " try:\n",
- " cwtmatr = cwt(h, ricker, [noise, ])\n",
- " except:\n",
- " return d\n",
- " marg = np.argmax(cwtmatr)\n",
- " pc = c[marg]\n",
- " # too large shift to be easily decernable via noise\n",
- " if pc > 10 or pc < -baseline_corr_noise_threshold:\n",
- " return d\n",
- " return d - pc\n",
- " \n",
- " def read_fun(filename, channel):\n",
- " \n",
- " infile = h5py.File(filename, \"r\", driver=\"core\")\n",
- " if rawversion == 2:\n",
- " count = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/count\".format(instrument, channel)])\n",
- " first = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/first\".format(instrument, channel)])\n",
- " last_index = int(first[count != 0][-1]+count[count != 0][-1])\n",
- " first_index = int(first[count != 0][0])\n",
- " else:\n",
- " status = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/status\".format(instrument, channel)])\n",
- " if np.count_nonzero(status != 0) == 0:\n",
- " return\n",
- " last = np.squeeze(infile[\"/INDEX/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/last\".format(instrument, channel)])\n",
- " last_index = int(last[status != 0][-1])\n",
- " first_index = int(last[status != 0][0])\n",
- " \n",
- " if limit_trains is not None:\n",
- " last_index = min(limit_trains*memory_cells+first_index, last_index)\n",
- " \n",
- " im = np.array(infile[\"/INSTRUMENT/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/data\".format(instrument, channel)][first_index:last_index,...]) \n",
- " carr = infile[\"/INSTRUMENT/{}_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/cellId\".format(instrument, channel)][first_index:last_index]\n",
- " cells = np.squeeze(np.array(carr))\n",
- " infile.close()\n",
- " \n",
- " \n",
- " if il_mode:\n",
- " ga = im[1::2, 0, ...]\n",
- " im = im[0::2, 0, ...].astype(np.float32)\n",
- " else:\n",
- " ga = im[:, 1, ...]\n",
- " im = im[:, 0, ...].astype(np.float32)\n",
- "\n",
- " im = np.rollaxis(im, 2)\n",
- " im = np.rollaxis(im, 2, 1)\n",
- "\n",
- " ga = np.rollaxis(ga, 2)\n",
- " ga = np.rollaxis(ga, 2, 1)\n",
- " return im, ga, cells\n",
- " \n",
- " offset_cor = xcal.OffsetCorrection(sensor_size, offset, nCells=memory_cells, blockSize=block_size, gains=[0,1,2])\n",
- " offset_cor.debug()\n",
- " \n",
- " hist_calc = xcal.HistogramCalculator(sensor_size, bins=2000, range=(0, 2000), blockSize=block_size)\n",
- " hist_calc.debug()\n",
- "\n",
- " hist_calc_uncorr = xcal.HistogramCalculator(sensor_size, bins=2000, range=(0, 2000), blockSize=block_size)\n",
- " hist_calc_uncorr.debug()\n",
- " \n",
- " \n",
- " for run in runs:\n",
- " for seq in sequences[:2]:\n",
- " \n",
- " fname = fbase.format(run, run.upper(), channel, seq)\n",
- " \n",
- " d, ga, c = read_fun(fname, channel)\n",
- " \n",
- " # we need to do proper gain thresholding\n",
- " vidx = (c < offset.shape[2]) & (c < thresholds.shape[2])\n",
- " c = c[vidx]\n",
- " d = d[:,:,vidx]\n",
- " ga = ga[:,:,vidx]\n",
- "\n",
- " g = np.zeros(ga.shape, np.uint8)\n",
- " g[...] = 2\n",
- " \n",
- " g[ga < thresholds[...,c,1]] = 1\n",
- " g[ga < thresholds[...,c,0]] = 0\n",
- " d = offset_cor.correct(d, cellTable=c, gainMap=g)\n",
- " for i in range(d.shape[2]):\n",
- " mn_noise = np.nanmean(noise[...,c[i]])\n",
- " d[...,i] = baseline_correct_via_noise(d[...,i],\n",
- " mn_noise,\n",
- " g[..., i])\n",
- " \n",
- " d *= np.moveaxis(pc[c,...], 0, 2)\n",
- " \n",
- " hist_calc_uncorr.fill(d)\n",
- " d = (d-gains[..., 1][...,None])/gains[..., 0][...,None]\n",
- " #d = d/gains[..., 0][...,None]\n",
- " hist_calc.fill(d) \n",
- " \n",
- " h, e, c, _ = hist_calc.get()\n",
- " hu = hist_calc_uncorr.get()\n",
- " return h, e, c, hu[0]\n",
- "\n",
- "inp = []\n",
- "idx = 0\n",
- "for i in modules:\n",
- " qm = \"Q{}M{}\".format(i//4+1, i%4+1)\n",
- " \n",
- " inp.append((i, offset_g[qm], thresholds_g[qm], res_gain[idx], noise_g[qm][...,0], pc_g[qm][0,...]))\n",
- " idx += 1\n",
- " \n",
- "p = partial(hist_single_module, fbase, runs, sequences, IL_MODE, profile, limit_trains, memory_cells, rawversion, instrument)\n",
- "#res = view.map_sync(p, inp)\n",
- "res = list(map(p, inp))\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "ExecuteTime": {
- "end_time": "2019-09-13T09:23:56.023378Z",
- "start_time": "2019-09-13T09:23:55.002226Z"
- }
- },
- "outputs": [],
- "source": [
- "from iminuit import Minuit\n",
- "from iminuit.util import make_func_code, describe\n",
- "from IPython.display import HTML, display\n",
- "import tabulate\n",
- "\n",
- "# fitting\n",
- "par_ests = {}\n",
- "par_ests[\"mu0\"] = 0\n",
- "par_ests[\"mu1\"] = 50\n",
- "par_ests[\"mu2\"] = 100\n",
- "par_ests[\"limit_mu0\"] = [-5, 5]\n",
- "par_ests[\"limit_mu1\"] = [35, 65]\n",
- "par_ests[\"limit_mu2\"] = [100, 150]\n",
- "par_ests[\"s0\"] = 10\n",
- "par_ests[\"s1\"] = 10\n",
- "par_ests[\"s2\"] = 10\n",
- "par_ests[\"limit_A0\"] = [1e5, 1e12]\n",
- "par_ests[\"limit_A1\"] = [1e4, 1e10]\n",
- "par_ests[\"limit_A2\"] = [1e3, 1e10]\n",
- "\n",
- "\n",
- "par_ests[\"throw_nan\"] = False\n",
- "par_ests[\"pedantic\"] = False\n",
- "par_ests[\"print_level\"] = 1\n",
- "\n",
- "def gaussian3(x, mu0, s0, A0, mu1, s1, A1, mu2, s2, A2):\n",
- " return (A0/np.sqrt(2*np.pi*s0**2)*np.exp(-0.5*((x-mu0)/s0)**2) +\n",
- " A1/np.sqrt(2*np.pi*s1**2)*np.exp(-0.5*((x-mu1)/s1)**2) +\n",
- " A2/np.sqrt(2*np.pi*s2**2)*np.exp(-0.5*((x-mu2)/s2)**2))\n",
- "\n",
- "\n",
- "f_sig = describe(gaussian3)[1:]\n",
- "\n",
- "class _Chi2Functor:\n",
- " def __init__(self, f, x, y, err):\n",
- " self.f = f\n",
- " self.x = x\n",
- " self.y = y\n",
- " self.err = err\n",
- " f_sig = describe(f)\n",
- " # this is how you fake function\n",
- " # signature dynamically\n",
- " self.func_code = make_func_code(\n",
- " f_sig[1:]) # docking off independent variable\n",
- " self.func_defaults = None # this keeps numpy.vectorize happy\n",
- "\n",
- " def __call__(self, *arg):\n",
- " # notice that it accept variable length\n",
- " # positional arguments\n",
- " # chi2 = sum((y-self.f(x,*arg))**2 for x,y in zip(self.x,self.y))\n",
- " return np.sum(((self.f(self.x, *arg) - self.y) ** 2) / self.err)\n",
- " \n",
- " \n",
- "d = []\n",
- "y_range_max = 0\n",
- "table = []\n",
- "headers = ['Module',\n",
- " 'FWHM (cor.) [ADU]', 'Separation (cor.) [$\\sigma$]',\n",
- " 'FWHM (uncor.) [ADU]', 'Separation (uncor.) [$\\sigma$]', \n",
- " 'Improvement'\n",
- " ]\n",
- "for i, r in enumerate(res):\n",
- " qm = \"Q{}M{}\".format(i//4+1, i%4+1)\n",
- " row = []\n",
- " row.append(qm)\n",
- " \n",
- " h, e, c, hu = r\n",
- " \n",
- " \n",
- " d.append({\n",
- " 'x': c, \n",
- " 'y': h,\n",
- " 'drawstyle': 'steps-mid',\n",
- " 'label': '{}: corrected'.format(qm),\n",
- " 'marker': '.',\n",
- " 'color': 'blue',\n",
- " })\n",
- " \n",
- " idx = (h > 0) & np.isfinite(h)\n",
- " h_non_zero = h[idx]\n",
- " c_non_zero = c[idx]\n",
- " par_ests[\"A0\"] = np.float(h[np.argmin(abs(c-0))])\n",
- " par_ests[\"A1\"] = np.float(h[np.argmin(abs(c-50))])\n",
- " par_ests[\"A2\"] = np.float(h[np.argmin(abs(c-100))])\n",
- " wrapped = _Chi2Functor(gaussian3, c_non_zero, h_non_zero,\n",
- " np.sqrt(h_non_zero))\n",
- " \n",
- " m = Minuit(wrapped, **par_ests)\n",
- " fmin = m.migrad()\n",
- " \n",
- " xt = np.arange(0, 200)\n",
- " \n",
- " yt = gaussian3(xt, m.values['mu0'], m.values['s0'], m.values['A0'],\n",
- " m.values['mu1'], m.values['s1'], m.values['A1'],\n",
- " m.values['mu2'], m.values['s2'], m.values['A2'])\n",
- "\n",
- " d.append({\n",
- " 'x': xt, \n",
- " 'y': yt,\n",
- " 'label': '{}: corrected (fit)'.format(qm),\n",
- " 'color': 'green',\n",
- " 'drawstyle': 'line',\n",
- " 'linewidth': 2\n",
- " })\n",
- " \n",
- " \n",
- " d.append({\n",
- " 'x': c, \n",
- " 'y': hu,\n",
- " 'drawstyle': 'steps-mid',\n",
- " 'label': '{}: uncorrected'.format(qm),\n",
- " 'marker': '.',\n",
- " 'color': 'grey',\n",
- " 'alpha': 0.5\n",
- " })\n",
- " \n",
- " row += [m.values['s1']*2.35, (m.values['mu1']-m.values['mu0'])/m.values['s1']]\n",
- " \n",
- " \n",
- " idx = (hu > 0) & np.isfinite(hu)\n",
- " h_non_zero = hu[idx]\n",
- " c_non_zero = c[idx]\n",
- " wrapped = _Chi2Functor(gaussian3, c_non_zero, h_non_zero,\n",
- " np.sqrt(h_non_zero))\n",
- " \n",
- " m = Minuit(wrapped, **par_ests)\n",
- " fmin = m.migrad()\n",
- " \n",
- " xt = np.arange(0, 200)\n",
- " \n",
- " yt = gaussian3(xt, m.values['mu0'], m.values['s0'], m.values['A0'],\n",
- " m.values['mu1'], m.values['s1'], m.values['A1'],\n",
- " m.values['mu2'], m.values['s2'], m.values['A2'])\n",
- "\n",
- " d.append({\n",
- " 'x': xt, \n",
- " 'y': yt,\n",
- " 'label': '{}: uncorrected (fit)'.format(qm),\n",
- " 'color': 'red',\n",
- " 'linewidth': 2\n",
- " })\n",
- " \n",
- " row += [m.values['s1']*2.35, (m.values['mu1']-m.values['mu0'])/m.values['s1']]\n",
- " \n",
- " row.append(\"{:0.2f} %\".format(100*(row[3]/row[1]-1)))\n",
- " \n",
- " y_range_max = max(y_range_max, np.max(h[c>25])*1.5)\n",
- " \n",
- " # output table\n",
- " table.append(row)\n",
- " \n",
- "fig = xana.simplePlot(d, y_log=False, figsize=\"2col\",\n",
- " aspect=2,\n",
- " x_range=(0, 200),\n",
- " legend='top-right-frame',\n",
- " y_range=(0, y_range_max),\n",
- " x_label='Energy (ADU)',\n",
- " y_label='Counts')\n",
- "\n",
- "display(HTML(tabulate.tabulate(table, tablefmt='html', headers=headers,\n",
- " numalign=\"right\", floatfmt=\"0.2f\")))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Python 3",
- "language": "python",
- "name": "python3"
- },
- "language_info": {
- "codemirror_mode": {
- "name": "ipython",
- "version": 3
- },
- "file_extension": ".py",
- "mimetype": "text/x-python",
- "name": "python",
- "nbconvert_exporter": "python",
- "pygments_lexer": "ipython3",
- "version": "3.6.7"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 1
-}
diff --git a/xfel_calibrate/notebooks.py b/xfel_calibrate/notebooks.py
index 54244492bb714f5d84ee8cbfcac949761e82aed7..7c73714ea8f52a6bdf061ae1edfe8417e1a7eaf7 100644
--- a/xfel_calibrate/notebooks.py
+++ b/xfel_calibrate/notebooks.py
@@ -21,8 +21,6 @@ notebooks = {
"FF": {
"notebook":
"notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_NBC.ipynb",
- "dep_notebooks": [
- "notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_Summary.ipynb"],
"concurrency": {"parameter": "modules",
"default concurrency": 16,
"cluster cores": 16},
@@ -172,7 +170,7 @@ notebooks = {
"notebook":
"notebooks/Jungfrau/Jungfrau_dark_analysis_all_gains_burst_mode_NBC.ipynb", # noqa
"concurrency": {"parameter": "karabo_da",
- "default concurrency": None,
+ "default concurrency": list(range(8)),
"cluster cores": 4},
},
"CORRECT": {