ImXPAD S540 detector at D2AM ============================ This tutorial corresponds to the calibration the goniometer an ImXPAD detector composed of 8 stripes of 7 modules, many of which are defective, on a goniometer. This detector is mounted on the goniometer 2theta arm at the D2AM beam-line, French CRG at the ESRF synchrotron. The raw data files are available at: http://www.silx.org/pub/pyFAI/gonio/D2AM-15/ .. code:: python %pylab nbagg .. parsed-literal:: Populating the interactive namespace from numpy and matplotlib .. code:: python import os, time start_time = time.time() import fabio, pyFAI from pyFAI.goniometer import GeometryTransformation, GoniometerRefinement, Goniometer from pyFAI.gui import jupyter .. code:: python #Download all data from silx.resources import ExternalResources #Specific for ESRF os.environ["http_proxy"] = "http://proxy.esrf.fr:3128" downloader = ExternalResources("pyFAI", "http://www.silx.org/pub/pyFAI/testimages", "PYFAI_DATA") all_files = downloader.getdir("LaB6_gonio_D2AM.tar.bz2") print("List of files downloaded:") for i in all_files: print(" "+os.path.basename(i)) detector_file = [i for i in all_files if i.endswith("D5Geom.h5")][0] images = [i for i in all_files if i.endswith(".edf")] npt_files = [i for i in all_files if i.endswith(".npt")] .. parsed-literal:: List of files downloaded: LaB6_gonio_D2AM 16Dec08D5_1763-rsz.edf 16Dec08D5_1730-rsz.npt 16Dec08D5_1725-rsz.npt 16Dec08D5_1749-rsz.npt 16Dec08D5_1735-rsz.npt 16Dec08D5_1763-rsz.npt 16Dec08D5_1729-rsz.npt 16Dec08D5_1728-rsz.npt 16Dec08D5_1791-rsz.npt D5Geom.h5 16Dec08D5_1725-rsz.edf 16Dec08D5_1742-rsz.edf 16Dec08D5_1726-rsz.npt 16Dec08D5_1784-rsz.npt 16Dec08D5_1777-rsz.npt 16Dec08D5_1791-rsz.edf 16Dec08D5_1729-rsz.edf 16Dec08D5_1770-rsz.npt 16Dec08D5_1742-rsz.npt 16Dec08D5_1770-rsz.edf 16Dec08D5_1726-rsz.edf 16Dec08D5_1784-rsz.edf 16Dec08D5_1749-rsz.edf 16Dec08D5_1730-rsz.edf 16Dec08D5_1728-rsz.edf 16Dec08D5_1727-rsz.npt 16Dec08D5_1756-rsz.npt 16Dec08D5_1756-rsz.edf 16Dec08D5_1735-rsz.edf 16Dec08D5_1727-rsz.edf 16Dec08D5_1777-rsz.edf .. code:: python #Definition of the detector and deplay of an image and its mask: d5 = pyFAI.detector_factory(detector_file) print("Detector shape: ",d5.shape) fimg = fabio.open(images[-1]) for k,v in fimg.header.items(): print(k, ": ", v) f, ax = subplots(1, 2) jupyter.display(d5.mask, label="mask", ax=ax[0]) jupyter.display(fimg.data, label=os.path.basename(fimg.filename), ax=ax[1]) .. parsed-literal:: Detector shape: (960, 578) EDF_DataBlockID : 0.Image.Psd EDF_BinarySize : 4439040 EDF_HeaderSize : 1536 ByteOrder : LowByteFirst DataType : DoubleValue Dim_1 : 578 Dim_2 : 960 Image : 0 HeaderID : EH:000000:000000:000000 Size : 4439040 VersionNumber : 1 Epoch : 1481330682.5452539921 det_sample_dist : 0 y_beam : 0 x_beam : 0 Lambda : 0.495938 offset : 0 count_time : 120 point_no : 52 scan_no : 906 preset : 0 col_end : 559 col_beg : 0 row_end : 959 row_beg : 0 counter_pos : 120 2733 166 15.3905 0 91.0146 15.3905 0 25 25 1777 7.45541e+09 0 0 0 174.822 0 0 15.3905 counter_mne : sec vct1 vct2 vct3 vct4 Imach pseudoC pfoil Emono Ecod img roi1 roi2 roi3 roi4 pico1 pico2 pico3 pico4 motor_pos : 52 0.078446 89.991 -89.9919 -0.0031 0.0021 57.1196 134.747 -32.9501 0.16656 -5 0.47558 -1.5 0 0 0 4.53604 0.1416 1.04 1.04022 1.04022 -4.4 -1.10211 -0.543725 -9.962 -14.038 -16.865 -7.195 24 -2.038 24.06 4.835 1 motor_mne : del eta chi phi nu mu keta kap kphi tsx tsy tsz rox roy tox toy mono inc1 courb courbb courbf omega khimono gamma su6 sd6 sf6 sb6 vg6 vo6 hg6 ho6 rien suffix : .edf prefix : 16Dec08D5_ dir : /users/opd02/raw run : 1777 title : CCD Image .. parsed-literal:: .. raw:: html .. parsed-literal:: .. code:: python # Define wavelength and create our "large" LaB6 calibrant wavelength = 0.495938 * 1e-10 from pyFAI.calibrant import get_calibrant LaB6 = get_calibrant("LaB6") LaB6.wavelength = wavelength print("2theta max: ", numpy.degrees(LaB6.get_2th()[-1])) print("Number of reflections: ", len(LaB6.get_2th())) .. parsed-literal:: 2theta max: 179.173497672 Number of reflections: 236 .. code:: python #Use a few manually calibrated images: npt_files.sort() print("Number of hand-calibrated images :",len(npt_files)) .. parsed-literal:: Number of hand-calibrated images : 15 .. code:: python # Definition of the goniometer translation function: # The detector rotates vertically, around the horizontal axis, i.e. rot2. # Rotation both around axis 1 and axis 2 are allowed goniotrans = GeometryTransformation(param_names = ["dist", "poni1", "poni2", "rot1", "rot2", "rot3", "scale1", "scale2" ], dist_expr="dist", poni1_expr="poni1", poni2_expr="poni2", rot1_expr="scale1 * pos +rot1", rot2_expr="scale2 * pos + rot2", rot3_expr="rot3") #Definition of the function reading the goniometer angle from the filename of the image. def get_angle(metadata): """Takes the angle from the first motor position and returns the angle of the goniometer arm""" return float(metadata["motor_pos"].split()[0]) print('filename', os.path.basename(fimg.filename), "angle:",get_angle(fimg.header)) .. parsed-literal:: filename 16Dec08D5_1777-rsz.edf angle: 52.0 .. code:: python # Definition of the geometry refinement: the parameter order is the same as the param_names rot3 = numpy.pi/2 scale1 = -numpy.pi/180 scale2 = 0 param = {"dist":0.5, "poni1":0.05, "poni2":0.05, "rot1":0, "rot2":0, "rot3": rot3, "scale1": scale1, "scale2": scale2, } #Defines the bounds for some variables bounds = {"dist": (0.2, 0.8), "poni1": (0, 0.1), "poni2": (0, 0.1), "rot1": (-0.1, 0.1), "rot2": (-0.1, 0.1), "rot3": (rot3, rot3), #strict bounds on rot3 #"scale1": (scale1, scale1), #"scale2": (scale2, scale2), } gonioref = GoniometerRefinement(param, #initial guess bounds=bounds, pos_function=get_angle, trans_function=goniotrans, detector=d5, wavelength=wavelength) print("Empty refinement object:", gonioref) #Let's populate the goniometer refinement object with all control point files: for fn in npt_files[:]: base = os.path.splitext(fn)[0] fimg = fabio.open(base + ".edf") basename = os.path.basename(base) sg =gonioref.new_geometry(basename, image=fimg.data, metadata=fimg.header, control_points=fn, calibrant=LaB6) print(basename, "Angle:", sg.get_position()) print("Filled refinement object:") print(gonioref) .. parsed-literal:: Empty refinement object: GoniometerRefinement with 0 geometries labeled: . 16Dec08D5_1725-rsz Angle: -0.003 16Dec08D5_1726-rsz Angle: 0.9998 16Dec08D5_1727-rsz Angle: 2.0 16Dec08D5_1728-rsz Angle: 2.9998 16Dec08D5_1729-rsz Angle: 4.0002 16Dec08D5_1730-rsz Angle: 4.9998 16Dec08D5_1735-rsz Angle: 10.0001 16Dec08D5_1742-rsz Angle: 16.9996 16Dec08D5_1749-rsz Angle: 24.0001 16Dec08D5_1756-rsz Angle: 30.9997 16Dec08D5_1763-rsz Angle: 37.9999 16Dec08D5_1770-rsz Angle: 44.9997 16Dec08D5_1777-rsz Angle: 52.0 16Dec08D5_1784-rsz Angle: 58.9995 16Dec08D5_1791-rsz Angle: 65.9999 Filled refinement object: GoniometerRefinement with 15 geometries labeled: 16Dec08D5_1725-rsz, 16Dec08D5_1726-rsz, 16Dec08D5_1727-rsz, 16Dec08D5_1728-rsz, 16Dec08D5_1729-rsz, 16Dec08D5_1730-rsz, 16Dec08D5_1735-rsz, 16Dec08D5_1742-rsz, 16Dec08D5_1749-rsz, 16Dec08D5_1756-rsz, 16Dec08D5_1763-rsz, 16Dec08D5_1770-rsz, 16Dec08D5_1777-rsz, 16Dec08D5_1784-rsz, 16Dec08D5_1791-rsz. .. code:: python # Initial refinement of the goniometer model with 5 dof gonioref.refine2() .. parsed-literal:: Cost function before refinement: 0.000307661216084 fun: 1.9855225499763415e-08 jac: array([ 2.40901503e-08, 3.04893444e-07, 4.13717251e-07, -1.75578498e-07, -9.57395104e-07, -4.24105195e-14, -5.24620326e-06, -7.15065733e-07]) message: 'Optimization terminated successfully.' nfev: 178 nit: 17 njev: 17 status: 0 success: True x: array([ 5.20464974e-01, 6.34779274e-02, 4.46077058e-02, 3.03313915e-03, 7.46909754e-03, 1.57079633e+00, -1.74719084e-02, -3.84262989e-04]) Cost function after refinement: 1.9855225499763415e-08 GonioParam(dist=0.52046497446866713, poni1=0.063477927428416264, poni2=0.04460770577505211, rot1=0.0030331391536378505, rot2=0.0074690975437223016, rot3=1.5707963267948966, scale1=-0.017471908444279767, scale2=-0.00038426298857720866) maxdelta on: dist (0) 0.5 --> 0.520464974469 .. parsed-literal:: array([ 5.20464974e-01, 6.34779274e-02, 4.46077058e-02, 3.03313915e-03, 7.46909754e-03, 1.57079633e+00, -1.74719084e-02, -3.84262989e-04]) .. code:: python width=3 height=int(ceil(len(gonioref.single_geometries) / width)) fig,ax = subplots(height, width,figsize=(10, 25)) for idx, sg in enumerate(gonioref.single_geometries.values()): sg.geometry_refinement.set_param(gonioref.get_ai(sg.get_position()).param) jupyter.display(sg=sg, ax=ax[idx//width, idx%width]) .. parsed-literal:: .. raw:: html .. parsed-literal:: /scisoft/users/jupyter/jupy34/lib/python3.4/site-packages/matplotlib/axes/_axes.py:545: UserWarning: No labelled objects found. Use label='...' kwarg on individual plots. warnings.warn("No labelled objects found. " .. code:: python # Final pass of refinement with all constrains removed, very fine refinement gonioref.bounds = None gonioref.refine2("slsqp", eps=1e-13, maxiter=10000, ftol=1e-12) .. parsed-literal:: Cost function before refinement: 1.98552254998e-08 fun: 1.9855093592191876e-08 jac: array([ 1.03563013e-08, 0.00000000e+00, 0.00000000e+00, 4.48001020e-08, -8.71252796e-07, 2.31610572e-09, -3.46915909e-05, -5.50968462e-07]) message: 'Optimization terminated successfully.' nfev: 15 nit: 1 njev: 1 status: 0 success: True x: array([ 5.20464974e-01, 6.34779274e-02, 4.46077058e-02, 3.03313914e-03, 7.46909773e-03, 1.57079633e+00, -1.74719008e-02, -3.84262868e-04]) Cost function after refinement: 1.9855093592191876e-08 GonioParam(dist=0.52046497446639939, poni1=0.063477927428416264, poni2=0.04460770577505211, rot1=0.0030331391438279576, rot2=0.0074690977345008092, rot3=1.5707963267943894, scale1=-0.01747190084785084, scale2=-0.00038426286793146559) maxdelta on: scale1 (6) -0.0174719084443 --> -0.0174719008479 .. parsed-literal:: array([ 5.20464974e-01, 6.34779274e-02, 4.46077058e-02, 3.03313914e-03, 7.46909773e-03, 1.57079633e+00, -1.74719008e-02, -3.84262868e-04]) .. code:: python #Create a MultiGeometry integrator from the refined geometry: angles = [] images = [] for sg in gonioref.single_geometries.values(): angles.append(sg.get_position()) images.append(sg.image) multigeo = gonioref.get_mg(angles) multigeo.radial_range=(0, 80) print(multigeo) .. parsed-literal:: MultiGeometry integrator with 15 geometries on (0, 80) radial range (2th_deg) and (-180, 180) azimuthal range (deg) .. code:: python # Integrate the whole set of images in a single run: res = multigeo.integrate1d(images, 10000) jupyter.plot1d(res) #Note the large number of peaks due to hot pixels .... .. parsed-literal:: .. raw:: html .. parsed-literal:: .. code:: python #Add hot pixels to the mask: pixel which are 15x more intense than the median in their ring. thres = 15 old_mask = d5.mask.astype("bool", copy=True) new_mask = d5.mask.astype("bool", copy=True) for ai,img in zip(multigeo.ais,images): b,a = ai.separate(img, 1000, restore_mask=0) b[old_mask] = 0 b[b<0] = 0 print(sum(b>thres*a)) new_mask = numpy.logical_or(new_mask, (b>thres*a)) print(" Size of old mask", sum(old_mask), " Size of new mask",sum(new_mask), " Number of pixel discarded", sum(new_mask)-sum(old_mask)) .. parsed-literal:: 4 4 4 4 4 4 4 4 4 4 4 3 3 3 3 Size of old mask 198678 Size of new mask 198683 Number of pixel discarded 5 .. code:: python # Update the mask for ai in multigeo.ais: ai.detector.mask = new_mask # Integrate the whole set of images in a single run: res2 = multigeo.integrate1d(images, 10000) ax = jupyter.plot1d(res, label="Before hot-pixel removal") ax.plot(*res2, label="After hot-pixel removal") ax.legend() .. parsed-literal:: .. raw:: html .. parsed-literal:: .. code:: python # Integrate the whole set of images in 2D: res2d = multigeo.integrate2d(images, 1000, 360) jupyter.plot2d(res2d) .. parsed-literal:: .. raw:: html .. parsed-literal:: .. code:: python print("Total execution time", time.time()-start_time) .. parsed-literal:: Total execution time 60.505094051361084