diff --git a/ImageD11/nbGui/S3DXRD/select_for_index_unknown.ipynb b/ImageD11/nbGui/S3DXRD/select_for_index_unknown.ipynb index efa07775..9e11d6cc 100644 --- a/ImageD11/nbGui/S3DXRD/select_for_index_unknown.ipynb +++ b/ImageD11/nbGui/S3DXRD/select_for_index_unknown.ipynb @@ -115,18 +115,9 @@ "outputs": [], "source": [ "# powder patterns\n", - "hsI, dsbinedges = np.histogram( cf_4d.ds, weights=cf_4d.sum_intensity, bins=np.linspace(0, cf_4d.ds.max(), 5000))\n", - "hs1, dsbinedges = np.histogram( cf_4d.ds, bins=np.linspace(0, cf_4d.ds.max(), 5000))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ + "pbins = 5000\n", + "hsI, dsbinedges = np.histogram( cf_4d.ds, weights=cf_4d.sum_intensity, bins=np.linspace(0, cf_4d.ds.max(), pbins))\n", + "hs1, dsbinedges = np.histogram( cf_4d.ds, bins=np.linspace(0, cf_4d.ds.max(), pbins))\n", "dsbincens = 0.5*(dsbinedges[1:]+dsbinedges[:-1])" ] }, @@ -166,7 +157,9 @@ "outputs": [], "source": [ "# These peaks are already indexed by the main phase\n", - "ringmask = ImageD11.peakselect.rings_mask( cf_4d, 0.005, 1.0, cell = ucell )" + "dstar_tol = 0.005 # range in dstar to mask\n", + "dsmax = 1.0 # max dstar (highest angle)\n", + "ringmask = ImageD11.peakselect.rings_mask( cf_4d, dstar_tol, 1.0, cell = ucell )" ] }, { @@ -178,7 +171,7 @@ "outputs": [], "source": [ "# cring = columnfile of things not indexed by the main rings\n", - "cring = cf_4d.copyrows( (cf_4d.ds < 1.)&( ~ringmask ))" + "cring = cf_4d.copyrows( (cf_4d.ds < dsmax)&( ~ringmask ))" ] }, { @@ -193,7 +186,7 @@ "ax.plot( dsbincens, hsI/hsI.max() , '-')\n", "ax.plot( ucell.ringds, np.full( len(ucell.ringds), 1e-5) , \"|\", ms=20, lw =1 )\n", "ax.plot( cring.ds, cring.sum_intensity/cring.sum_intensity.max(), \".\", ms=1 )\n", - "ax.set(yscale='log', xlim=(0.1,1), ylim=(1e-5,1))" + "ax.set(yscale='log', xlim=(0.1,1), ylim=(1e-5,1));" ] }, { @@ -240,6 +233,20 @@ }, "outputs": [], "source": [ + "# Position where the beam goes through the rotation axis\n", + "y0 = 0.0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Interactive plot to select a peak\n", + "\n", "fig, ax = plt.subplots(1,2)\n", "ax[0].scatter( cpk.omega, cpk.dty, \n", " c=np.log(cpk.sum_intensity) * np.sign( cpk.yl ), \n", @@ -247,7 +254,6 @@ "\n", "selected, = ax[0].plot( [], [], \"o\", mfc = 'none')\n", "\n", - "\n", "ax[1].scatter( cring.omega, cring.dty,\n", " c=np.log(cring.sum_intensity) * np.sign( cring.yl ), \n", " s=4, cmap='RdBu')\n", @@ -255,7 +261,8 @@ "\n", "def find_pair( idx, cf ):\n", " \"\"\"\n", - " Given a peak idx, locate the Friedel pair in cf\n", + " Given a peak idx, locate the Friedel pair in colfile\n", + " Matches on gvector only (not intensity, etc)\n", " \"\"\"\n", " dg2 = (cf.gx[idx] + cf.gx)**2 + (cf.gy[idx] + cf.gy)**2 + (cf.gz[idx] + cf.gz)**2 \n", " pair = np.argmin( dg2 )\n", @@ -264,32 +271,42 @@ "om = np.linspace( cring.omega.min(), cring.omega.max(), 90 )\n", "fitline, = ax[1].plot( om, np.zeros_like(om), '-')\n", "\n", - "### GLOBALS\n", - "p1 = None\n", - "p2 = None\n", - "xy = None\n", - "peak_ycalc = None\n", + "### GLOBALS, also y0\n", + "p1 = None # index of first peak in cring\n", + "p2 = None # index of second peak in cring\n", + "xy = None # fitted position\n", + "peak_ycalc = None # dty computed\n", "\n", "def fit_pair( cf, i, j ):\n", " # fits a sine wave\n", - " global xy, p1, p2, peak_ycalc\n", + " global xy, p1, p2, peak_ycalc, y0\n", " p1 = i\n", " p2 = j\n", + " # sin/cos omega\n", " ci = np.cos( np.radians( cf.omega[i] ) )\n", " cj = np.cos( np.radians( cf.omega[j] ) )\n", " si = np.sin( np.radians( cf.omega[i] ) )\n", " sj = np.sin( np.radians( cf.omega[j] ) )\n", - " yi = cf.dty[i]\n", - " yj = cf.dty[j]\n", - " # yi = xp.ci + yp.si\n", - " # yj = xp.cj + yp.sj\n", " #\n", + " # ImageD11.sinograms.geometry:\n", + " # def dtycalc_sincos(sinomega, cosomega, x, y, y0):\n", + " # ....\n", + " # return y0 - x * sinomega - y * cosomega\n", + " #\n", + " # yi = y0 - xp.si - yp.ci\n", + " # yj = y0 - xp.sj - yp.cj\n", + " #\n", + " yi = cf.dty[i] - y0\n", + " yj = cf.dty[j] - y0\n", + " #\n", + " # Only two peaks, cannot fit x/y/y0\n", " # Y = R.x\n", - " R = np.array( [ [ ci, si ], [cj, sj] ] )\n", - " xy = np.linalg.inv(R).dot( [ yi, yj ] )\n", - " ycalc = np.cos(np.radians(om))*xy[0]+np.sin(np.radians(om))*xy[1]\n", + " R = np.array( [ [ -si, -ci ], \n", + " [ -sj, -cj] ] )\n", + " xy = x, y = np.linalg.inv(R).dot( [ yi, yj ] )\n", + " ycalc = ImageD11.sinograms.geometry.x_y_y0_omega_to_dty(om, x, y, y0)\n", " fitline.set_ydata( ycalc )\n", - " peak_ycalc = np.cos(np.radians(cring.omega))*xy[0]+np.sin(np.radians(cring.omega))*xy[1]\n", + " peak_ycalc = ImageD11.sinograms.geometry.x_y_y0_omega_to_dty(cring.omega, x, y, y0)\n", " \n", "def onpick( evt ):\n", " art = evt.artist\n", @@ -300,6 +317,7 @@ " selected.set_ydata( [cpk.dty[ind[-1]],cpk.dty[pair]] )\n", " ax[0].set( title = f'{ind} {pair}' )\n", " fig.canvas.draw_idle()\n", + "\n", "ax[1].set(title=\"Pick a peak on the other plot\")\n", "fig.canvas.mpl_connect( 'pick_event', onpick );" ] @@ -312,12 +330,13 @@ }, "outputs": [], "source": [ + "dtyrange = 5.\n", "print(\"Selected points are\",p1,p2)\n", "print(\"Omega, dty\")\n", "print(cpk.omega[p1],cpk.dty[p1])\n", "print(cpk.omega[p2],cpk.dty[p2])\n", "f,ax = plt.subplots()\n", - "ax.hist( cring.dty - peak_ycalc, bins = np.arange(-2,2,ds.ystep) );\n", + "ax.hist( cring.dty - peak_ycalc, bins = np.arange(-dtyrange,dtyrange,ds.ystep) );\n", "ax.set(title='dty error histogram', xlabel='dty', ylabel='frequency');" ] }, @@ -330,8 +349,9 @@ "outputs": [], "source": [ "# now we select the peaks from the sinogram based on that position in space\n", + "dty_cut = 0.4 # from the last plot\n", "\n", - "ymask = abs(peak_ycalc - cring.dty)<0.4\n", + "ymask = abs(peak_ycalc - cring.dty)=0:\n", + " \n", + " serial = silx.io.get_data( f\"silx:{masterfile}::{darkscan}/instrument/{detector}/camera_settings/camera_serial\" )\n", + " print(\"Detector serial\", serial )\n", + " if serial == \"21\":\n", + " flat = silx.io.get_data( f\"fabio:/data/id11/3dxrd/inhouse/Frelon21/F21_flat_july18.edf\" )\n", + " e2dxfile = \"/data/id11/3dxrd/inhouse/Frelon21/F21dxnew.edf\"\n", + " e2dyfile = \"/data/id11/3dxrd/inhouse/Frelon21/F21dynew.edf\"\n", + " pixelsize = float(fabio.open(e2dxfile).header['PSize_1'].split()[0])*1e6 # microns\n", + " print('Pixel size', pixelsize )\n", + " if serial == \"36\": \n", + " flat = silx.io.get_data( f\"fabio:/data/id11/3dxrd/inhouse/Frelon36/F36_Nov2023.edf\" )\n", + " flat = flat / flat[500:-500, 500:-500].mean()\n", + " e2dxfile = \"fabio:/data/id11/3dxrd/inhouse/Frelon36/f36dx.edf\"\n", + " e2dyfile = f\"fabio:/data/id11/3dxrd/inhouse/Frelon36/f36dy.edf\"\n", + " pixelsize = float(fabio.open(e2dxfile).header['PSize_1'].split()[0])*1e6 # microns\n", + " if serial == \"29\":\n", + " flat = silx.io.get_data( f\"fabio:/data/id11/3dxrd/inhouse/Frelon4M/c29_Agfilt_flat-dark.edf\" )\n", + " e2dxfile = \"fabio:/data/id11/3dxrd/inhouse/Frelon4M/F4M_EO_dx.edf\"\n", + " e2dyfile = \"fabio:/data/id11/3dxrd/inhouse/Frelon4M/F4M_EO_dy.edf\"\n", + " pixelsize = float(fabio.open(e2dxfile).header['PSize_1'].split()[0])*1e6 # microns\n", + " \n", + " print( \"Average flat, should be close to 1 : \",flat[100:-100, 100:-100].mean() )\n", + " \n", + " cor = ( frames - dark[ None, :, :] ).mean( axis = 0, dtype = np.float32 ) / flat\n", + " \n", + "else: # eiger\n", + " \n", + " cor = frames.mean( axis = 0, dtype = np.float32 )\n", + " \n", + " e2dxfile = \"/data/id11/nanoscope/Eiger/e2dx_E-08-0144_20240205.edf\"\n", + " e2dyfile = \"/data/id11/nanoscope/Eiger/e2dy_E-08-0144_20240205.edf\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2d4d3e39-6615-41da-9a49-7c3687306957", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def guess_median_bg( img, size=91 ):\n", + " m0 = scipy.ndimage.median_filter( cor, (size,1))\n", + " m1 = scipy.ndimage.median_filter( cor, (1,size))\n", + " bg = np.where(m0 lowcut )\n", + "\n", + "bg = guess_median_bg( cor )\n", + "\n", + "signal = ((cor - bg )*mask).clip(0,None)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "631b7295-1b53-4fac-994d-62967e242662", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "f,(a,b) = plt.subplots( 2, 2, figsize=(8, 6), constrained_layout=True )\n", + "f.colorbar( a[1].imshow( cor * mask ) )\n", + "f.colorbar( b[1].imshow( signal ) )\n", + "a[0].plot( cor[ 1024, : ] )\n", + "a[0].plot( cor[ :, 1024 ] )\n", + "a[0].set(xlabel='pixel', ylabel='intensity')\n", + "b[0].plot( signal[ 1024, : ] )\n", + "b[0].plot( signal[ :, 1024 ] )\n", + "b[0].set(xlabel='pixel', ylabel='intensity', ylim=(0, None) );" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0dd443d3-9be7-4c84-a1cc-502d3c4aab0a", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Old way, using command line:\n", + "# !powderimagetopeaks.py CeO2_for_S5_blue_rot_clean.edf CeO2_for_S5_blue_rot_clean0000.edf CeO2_for_S5_blue_rot_clean0001.edf 1024 1024\n", + "def powderimagetopeaks( img, ci = 1024, cj = 1024 ):\n", + " # cuts into 1 degree bins\n", + " si, sj = img.shape\n", + " i, j = np.mgrid[0:si,0:sj]\n", + " phi = np.arctan2( i - ci, j - cj ) * 361 / np.pi\n", + " return (phi%2).astype(int) == 0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dbc62cc4-03aa-4f33-a4a8-431b4e6fa00d", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# old way:\n", + "# !peaksearch.py -n CeO2_for_S5_blue_rot_clean -f 0 -l 1 -p N -s /data/id11/3dxrd/inhouse/Frelon36/frelon36.spline -t 5000 -t 10000 -t 20000 -t 40000 -o CeO2_for_S5_blue_rot_clean.spt " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ebe6849c-0b51-4ec2-95b8-2caefc9b75c5", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "worker = ImageD11.frelon_peaksearch.worker(None, None)\n", + "pkmsk = powderimagetopeaks( signal )\n", + "pks = np.concatenate( (\n", + " worker.peaksearch( signal * pkmsk ), # even degrees\n", + " worker.peaksearch( signal * (1-pkmsk)), # odd degrees\n", + " ) , axis = 0 )\n", + "pks.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d7cefae-e1f3-4bad-9a07-454b592a8ccc", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from ImageD11.cImageD11 import s_1, s_I, s_sI, s_fI\n", + "colfile = ImageD11.columnfile.colfile_from_dict( {\n", + " \"s_raw\" : pks[ :, s_sI ] / pks[ :, s_I ],\n", + " \"f_raw\" : pks[ :, s_fI ] / pks[ :, s_I ],\n", + " \"Number_of_pixels\" : pks[ :, s_1 ], \n", + " \"sum_intensity\" : pks[ :, s_I ],\n", + " \"omega\" : np.zeros( len(pks), float ),\n", + "} )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1fa38247-7d08-45cd-aeec-4618ac178917", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "f, a = plt.subplots(1,2, figsize=(8,5), constrained_layout=True)\n", + "a[0].plot( colfile.f_raw, colfile.s_raw, \".\", ms = 1)\n", + "a[0].set( ylabel=\"detector slow direction\", xlabel=\"detector fast direction\", aspect='equal')\n", + "a[1].plot( colfile.Number_of_pixels, colfile.sum_intensity, \".\", ms =1 )\n", + "a[1].set( xlabel=\"Number of pixels\", ylabel=\"sum intensity\", yscale='log', xscale='log');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2a543865-9505-4b25-bab3-f130546bcbff", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# old way\n", + "# !spatialfix.py frelon36_spline_dx.edf frelon36_spline_dy.edf CeO2_for_S5_blue_rot_clean_t5000.flt CeO2_for_S5_blue_rot_clean_spline_t5000.flt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "abec9195-69d1-44f8-9238-e549ee68fbf6", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "colfile = ImageD11.blobcorrector.eiger_spatial( e2dxfile, e2dyfile )(colfile)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bad00d35-5618-4a79-969b-ef2328e74447", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# input data from peaksearch and starting parameter file:\n", + "trans=ImageD11.transformer.transformer()\n", + "trans.colfile = colfile\n", + "trans.setxyomcols(\"sc\", \"fc\", \"omega\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "62bc178a-d49a-4242-b4a6-9566c324fd6a", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "distance = 210 * 1e3 # micron\n", + "wavelength = 12.3984 / energy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ce7a07c4-b353-4be7-989f-22a8150848f9", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "wavelength, energy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f9ba5882-5236-47bf-b85c-0982c4db89c0", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "trans.parameterobj.set( 'y_size', pixelsize )\n", + "trans.parameterobj.set( 'z_size', pixelsize )\n", + "trans.parameterobj.set( 'distance', distance )\n", + "trans.parameterobj.set( 'wavelength', wavelength )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "639cd3aa-b35d-41a7-acb5-335042c68fc0", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# unit cell\n", + "a = 5.4115260\n", + "[ trans.parameterobj.set( f'cell__{x}', a ) for x in 'abc']\n", + "[ trans.parameterobj.set( f'cell_{x}', 90 ) for x in ('alpha','beta','gamma')]\n", + "trans.parameterobj.set('cell_lattice_[P,A,B,C,I,F,R]', \"F\" )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20a04c19-5065-4edd-8241-df851a686c6d", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def parCallBack( arg ):\n", + " trans.parameterobj.parameters.update( { arg['owner'].description : arg['new'] } )\n", + " drawPlot()\n", + " \n", + "def fixVaryCallBack( arg ):\n", + " name = arg['owner'].description.split(\" \")[1]\n", + " vars = trans.getvars()\n", + " if arg.new and not (name in vars):\n", + " vars.append( name )\n", + " if name in vars and not arg.new:\n", + " vars.remove(name)\n", + " trans.parameterobj.set_varylist(vars)\n", + " \n", + "def fitCallBack(arg):\n", + " \"\"\" fit call back - runs fit \"\"\"\n", + " lo, hi = ax1.get_xlim()\n", + " trans.fit( lo, hi )\n", + " for i, pname in enumerate(vars):\n", + " layout[i,0].value = trans.parameterobj.get(pname)\n", + " drawPlot()\n", + " \n", + "def drawPlot():\n", + " tth, eta = trans.compute_tth_eta()\n", + " pt1.set_data( tth, eta )\n", + " fig.canvas.update\n", + "\n", + "# Things to be edited in the UI\n", + "vars = \"y_center z_center distance tilt_y tilt_z\".split()\n", + "steps = (1, 1, 100, 0.01, 0.01, 0)\n", + "nv = len(vars)\n", + "\n", + "# Draw the widgets:\n", + "layout = ipywidgets.GridspecLayout(nv+1,2)\n", + "for i,( pname, pstep ) in enumerate( zip( vars, steps ) ) :\n", + " layout[i,0] = ipywidgets.FloatText( description=pname, \n", + " value = trans.parameterobj.parameters.get(pname),\n", + " step=pstep)\n", + " layout[i,0].observe( parCallBack , names='value' )\n", + " layout[i,1] = ipywidgets.ToggleButton( description=\"Vary \"+pname, \n", + " value = pname in trans.getvars() )\n", + " layout[i,1].observe( fixVaryCallBack, names='value' )\n", + " \n", + "layout[nv,0] = ipywidgets.FloatText( description='fit_tolerance', \n", + " value = trans.parameterobj.parameters.get(\"fit_tolerance\"), step=0,)\n", + "layout[nv,0].observe( parCallBack , names='value' )\n", + "\n", + "layout[nv,1] = ipywidgets.Button(description=\"Run Fit (blocks)\")\n", + "layout[nv,1].on_click( fitCallBack )\n", + "\n", + "\n", + "# Draw plot\n", + "fig1 = plt.figure(figsize=(9,6))\n", + "ax1 = fig1.add_subplot()\n", + "tth, eta = trans.compute_tth_eta()\n", + "trans.addcellpeaks()\n", + "pt1, = ax1.plot( tth, eta, \",\")\n", + "ax1.set(xlabel=\"tth\", ylabel=\"eta\")\n", + "ax1.plot( trans.theorytth, [0,]*len(trans.theorytth), \"r|\", ms=360, alpha=0.2 )\n", + "# Add controls\n", + "display(layout)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36bbea15-894b-4424-abe9-1e34702da583", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import pprint\n", + "pprint.pprint( trans.parameterobj.parameters )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c0c4bd4e-b7d6-4234-8442-283f314c6a54", + "metadata": {}, + "outputs": [], + "source": [ + "trans.parameterobj.saveparameters('CeO2_fitted.par')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f2a17baf-133e-4a87-b900-8448bccc0eee", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "12.3985 / trans.parameterobj.get('wavelength')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a29f6ee3-95f5-4e97-b03b-4fdd6deb4c71", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (main)", + "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.11.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}