From 0587020a1091f86caabf67c4ab519a06da15347b Mon Sep 17 00:00:00 2001
From: Egor Sobolev <egor.sobolev@xfel.eu>
Date: Mon, 19 Aug 2024 21:04:40 +0200
Subject: [PATCH] Initial implementation of center refinement using powder
diffraction data
---
nb/powder.ipynb | 164 +++
pyproject.toml | 2 +
src/geomtools/powder/__init__.py | 2 +
src/geomtools/powder/calibrants/AgBh.D | 52 +
src/geomtools/powder/calibrants/LaB6.D | 1444 +++++++++++++++++++
src/geomtools/powder/calibrants/__init__.py | 16 +
src/geomtools/powder/find_center.py | 124 ++
src/geomtools/powder/image_tools.py | 100 ++
src/geomtools/powder/misc.py | 69 +
src/geomtools/powder/powder.py | 146 ++
src/geomtools/powder/powdersum.py | 356 +++++
src/geomtools/powder/refine.py | 75 +
src/geomtools/powder/shmemarray.py | 144 ++
13 files changed, 2694 insertions(+)
create mode 100644 nb/powder.ipynb
create mode 100644 src/geomtools/powder/__init__.py
create mode 100644 src/geomtools/powder/calibrants/AgBh.D
create mode 100644 src/geomtools/powder/calibrants/LaB6.D
create mode 100644 src/geomtools/powder/calibrants/__init__.py
create mode 100644 src/geomtools/powder/find_center.py
create mode 100644 src/geomtools/powder/image_tools.py
create mode 100644 src/geomtools/powder/misc.py
create mode 100644 src/geomtools/powder/powder.py
create mode 100644 src/geomtools/powder/powdersum.py
create mode 100644 src/geomtools/powder/refine.py
create mode 100644 src/geomtools/powder/shmemarray.py
diff --git a/nb/powder.ipynb b/nb/powder.ipynb
new file mode 100644
index 0000000..82bd3bd
--- /dev/null
+++ b/nb/powder.ipynb
@@ -0,0 +1,164 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c134ff02-0eba-4146-803b-2b69f9ccaef9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "runno = 189\n",
+ "propno = 900463\n",
+ "\n",
+ "ref_geom_fn = \"/gpfs/exfel/exp/SPB/202430/p900429/usr/Shared/geom/p5681_r64_67_refined.geom\" # initial geometry file name"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "0bdf665e-9ad7-4a40-a59e-23b428312f4f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# !geomtools-powder-center -g {ref_geom_fn} -c AgBh ../../../nb/powder_sum_p{propno:06d}_r{runno:04d}.h5"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "504eca33-7f3d-4b7e-9ce1-9edcadff6387",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import os\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "import h5py\n",
+ "\n",
+ "from extra_data import open_run\n",
+ "from extra_geom import agipd_asic_seams, AGIPD_1MGeometry\n",
+ "from extra_geom.motors import AGIPD_1MMotors\n",
+ "from extra.components import AGIPD1MQuadrantMotors\n",
+ "\n",
+ "from geomtools.powder.calibrants import get_calibrant\n",
+ "from geomtools.powder.misc import agipd_borders\n",
+ "from geomtools.powder import PowderDiffraction"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "33f4f527-d61e-45c1-8301-1be7c6aa8edf",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fn = f\"powder_sum_p{propno:06d}_r{runno:04d}.h5\"\n",
+ "if not os.path.exists(fn):\n",
+ " !srun -p exfel geomtools-powder-sum -p {propno} -r {runno} -D ^scratch/pycal/proc -m ^usr/Shared/mask/mask_hvoff_20240807.h5 -l agipd1m SPB_DET_AGIPD1M-1"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e9c147f0-e011-46cb-bf0f-e0ef01a6d8a4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "with h5py.File(fn, \"r\") as f:\n",
+ " # conditions\n",
+ " clen = f[\"powderSum/conditions/cameraLen\"][()]\n",
+ " photon_en = f[\"powderSum/conditions/photonEnergy\"][()]\n",
+ " lmd = f[\"powderSum/conditions/waveLength\"][()]\n",
+ " detector_id = f[\"powderSum/conditions/detectorId\"][()].decode()\n",
+ "\n",
+ " # average image\n",
+ " img_mean = f[\"powderSum/image/mean\"][:]\n",
+ " img_std = f[\"powderSum/image/std\"][:]\n",
+ " img_count = f[\"powderSum/image/count\"][:]\n",
+ " num_frames = f[\"powderSum/image/numFrames\"][:]\n",
+ " mask = f[\"powderSum/image/mask\"][:].astype(bool)\n",
+ " modules = f[\"powderSum/image/modules\"][:]\n",
+ "\n",
+ "print(\"Photon energy\", photon_en, \"(keV)\")\n",
+ "print(\"Wave length\", lmd, \"(nm)\")\n",
+ "print(\"Camera length\", clen, \"(m)\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8fdf43d7-db13-4493-a030-3132337e1512",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# geometry\n",
+ "run = open_run(propno, runno, data=\"raw\")\n",
+ "ref_geom = AGIPD_1MGeometry.from_crystfel_geom(ref_geom_fn)\n",
+ "tracker = AGIPD_1MMotors(ref_geom)\n",
+ "motors = AGIPD1MQuadrantMotors(run)\n",
+ "geom = tracker.geom_at(motors.most_frequent_positions())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "332a77f6-b3b7-4e13-b6e4-5838391949df",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "img = np.copy(img_mean)\n",
+ "geom.plot_data(img, norm=\"log\", interpolation=\"none\")\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "83dd3cb1-00d5-4337-80c5-e83d9eeec0b1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "rings = get_calibrant(\"AgBh\")\n",
+ "borders = agipd_borders()\n",
+ "pw = PowderDiffraction(geom, rings, clen, lmd,\n",
+ " border_mask=~borders, make_shadow_mask=True)\n",
+ "pw.fit(img_mean, mask, num_frames)\n",
+ "pw.refinement_info()\n",
+ "pw.draw()\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1bd08493-2562-4f92-82aa-eb866d5af791",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pw.draw_radial_profile(img_mean, mask)\n",
+ "plt.show()"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "egor-dev",
+ "language": "python",
+ "name": "egor-dev"
+ },
+ "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.9.16"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/pyproject.toml b/pyproject.toml
index 60b9735..b3af3ba 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -31,3 +31,5 @@ dependencies = [
[project.scripts]
geomtools-push = "geomtools.sfx.report:push_geometry"
geomtools-refine-sfx = "geomtools.sfx.optimiser:refine"
+geomtools-powder-sum = "geomtools.powder.powdersum:main"
+geomtools-powder-center = "geomtools.powder.find_center:main"
diff --git a/src/geomtools/powder/__init__.py b/src/geomtools/powder/__init__.py
new file mode 100644
index 0000000..5cfdbda
--- /dev/null
+++ b/src/geomtools/powder/__init__.py
@@ -0,0 +1,2 @@
+from .image_tools import make_shadow_mask # noqa E401
+from .powder import PowderDiffraction # noqa E401
diff --git a/src/geomtools/powder/calibrants/AgBh.D b/src/geomtools/powder/calibrants/AgBh.D
new file mode 100644
index 0000000..48cc5a2
--- /dev/null
+++ b/src/geomtools/powder/calibrants/AgBh.D
@@ -0,0 +1,52 @@
+# Calibrant: Silver Behenate (AgBh)
+# Pseudocrystal a=inf b=inf c=58.380
+# Ref: doi:10.1107/S0021889899012388
+5.83800000e+01 # (0,0,1)
+2.91900000e+01 # (0,0,2)
+1.94600000e+01 # (0,0,3)
+1.45950000e+01 # (0,0,4)
+1.16760000e+01 # (0,0,5)
+9.73000000e+00 # (0,0,6)
+8.34000000e+00 # (0,0,7)
+7.29750000e+00 # (0,0,8)
+6.48666667e+00 # (0,0,9)
+5.83800000e+00 # (0,0,10)
+5.30727273e+00 # (0,0,11)
+4.86500000e+00 # (0,0,12)
+4.49076923e+00 # (0,0,13)
+4.17000000e+00 # (0,0,14)
+3.89200000e+00 # (0,0,15)
+3.64875000e+00 # (0,0,16)
+3.43411765e+00 # (0,0,17)
+3.24333333e+00 # (0,0,18)
+3.07263158e+00 # (0,0,19)
+2.91900000e+00 # (0,0,20)
+2.78000000e+00 # (0,0,21)
+2.65363636e+00 # (0,0,22)
+2.53826087e+00 # (0,0,23)
+2.43250000e+00 # (0,0,24)
+2.33520000e+00 # (0,0,25)
+2.24538462e+00 # (0,0,26)
+2.16222222e+00 # (0,0,27)
+2.08500000e+00 # (0,0,28)
+2.01310345e+00 # (0,0,29)
+1.94600000e+00 # (0,0,30)
+1.88322581e+00 # (0,0,31)
+1.82437500e+00 # (0,0,32)
+1.76909091e+00 # (0,0,33)
+1.71705882e+00 # (0,0,34)
+1.66800000e+00 # (0,0,35)
+1.62166667e+00 # (0,0,36)
+1.57783784e+00 # (0,0,37)
+1.53631579e+00 # (0,0,38)
+1.49692308e+00 # (0,0,39)
+1.45950000e+00 # (0,0,40)
+1.42390244e+00 # (0,0,41)
+1.39000000e+00 # (0,0,42)
+1.35767442e+00 # (0,0,43)
+1.32681818e+00 # (0,0,44)
+1.29733333e+00 # (0,0,45)
+1.26913043e+00 # (0,0,46)
+1.24212766e+00 # (0,0,47)
+1.21625000e+00 # (0,0,48)
+1.19142857e+00 # (0,0,49)
diff --git a/src/geomtools/powder/calibrants/LaB6.D b/src/geomtools/powder/calibrants/LaB6.D
new file mode 100644
index 0000000..ffaf5ba
--- /dev/null
+++ b/src/geomtools/powder/calibrants/LaB6.D
@@ -0,0 +1,1444 @@
+# Calibrant: Lanthanum hexaboride (LaB6)
+# Primitive cubic cell a=4.1568 b=4.1568 c=4.1568 alpha=90.000 beta=90.000 gamma=90.000
+# Ref: NIST standard reference material 660c
+4.15682600 # (1, 0, 0) 6
+2.93931985 # (1, 1, 0) 12
+2.39994461 # (1, 1, 1) 8
+2.07841300 # (2, 0, 0) 6
+1.85898910 # (2, 1, 0) 24
+1.69701711 # (2, 1, 1) 24
+1.46965993 # (2, 2, 0) 12
+1.38560867 # (3, 0, 0) 30
+1.31450380 # (3, 1, 0) 24
+1.25333020 # (3, 1, 1) 24
+1.19997231 # (2, 2, 2) 8
+1.15289610 # (3, 2, 0) 24
+1.11095848 # (3, 2, 1) 48
+1.03920650 # (4, 0, 0) 6
+1.00817839 # (4, 1, 0) 48
+0.97977328 # (4, 1, 1) 36
+0.95364129 # (3, 3, 1) 24
+0.92949455 # (4, 2, 0) 24
+0.90709380 # (4, 2, 1) 48
+0.88623828 # (3, 3, 2) 24
+0.84850855 # (4, 2, 2) 24
+0.83136520 # (5, 0, 0) 30
+0.81522065 # (5, 1, 0) 72
+0.79998154 # (5, 1, 1) 32
+0.77190321 # (5, 2, 0) 72
+0.75892912 # (5, 2, 1) 48
+0.73482996 # (4, 4, 0) 12
+0.72361053 # (5, 2, 2) 48
+0.71288978 # (5, 3, 0) 48
+0.70263184 # (5, 3, 1) 48
+0.69280433 # (6, 0, 0) 30
+0.68337798 # (6, 1, 0) 24
+0.67432622 # (6, 1, 1) 72
+0.65725190 # (6, 2, 0) 24
+0.64918715 # (6, 2, 1) 96
+0.64141218 # (5, 4, 1) 48
+0.63391002 # (5, 3, 3) 24
+0.62666510 # (6, 2, 2) 24
+0.61966303 # (6, 3, 0) 72
+0.61289056 # (6, 3, 1) 48
+0.59998615 # (4, 4, 4) 8
+0.59383229 # (7, 0, 0) 54
+0.58786397 # (7, 1, 0) 84
+0.58207207 # (7, 1, 1) 48
+0.57644805 # (6, 4, 0) 24
+0.57098396 # (7, 2, 0) 72
+0.56567237 # (7, 2, 1) 96
+0.55547924 # (6, 4, 2) 48
+0.55058505 # (7, 2, 2) 48
+0.54581799 # (7, 3, 0) 24
+0.54117265 # (7, 3, 1) 72
+0.53222703 # (6, 5, 0) 72
+0.52791743 # (7, 3, 2) 96
+0.51960325 # (8, 0, 0) 6
+0.51559081 # (8, 1, 0) 96
+0.51166991 # (8, 1, 1) 96
+0.50783712 # (7, 3, 3) 24
+0.50408920 # (8, 2, 0) 48
+0.50042304 # (8, 2, 1) 96
+0.49683574 # (6, 5, 3) 48
+0.48988664 # (8, 2, 2) 36
+0.48651968 # (8, 3, 0) 48
+0.48322121 # (8, 3, 1) 120
+0.47998892 # (7, 5, 1) 56
+0.47682064 # (6, 6, 2) 24
+0.47371429 # (8, 3, 2) 96
+0.47066786 # (7, 5, 2) 48
+0.46474728 # (8, 4, 0) 24
+0.46186956 # (9, 0, 0) 102
+0.45904464 # (9, 1, 0) 48
+0.45627093 # (9, 1, 1) 72
+0.45354690 # (8, 4, 2) 48
+0.45087108 # (9, 2, 0) 48
+0.44824208 # (9, 2, 1) 120
+0.44311914 # (6, 6, 4) 24
+0.44062267 # (9, 2, 2) 144
+0.43816793 # (9, 3, 0) 120
+0.43575377 # (9, 3, 1) 48
+0.43104278 # (8, 5, 2) 48
+0.42874387 # (9, 3, 2) 96
+0.42425428 # (8, 4, 4) 24
+0.42206173 # (9, 4, 0) 48
+0.41990284 # (9, 4, 1) 108
+0.41777673 # (9, 3, 3) 72
+0.41568260 # (10, 0, 0) 30
+0.41361965 # (10, 1, 0) 168
+0.41158711 # (10, 1, 1) 48
+0.40761032 # (10, 2, 0) 72
+0.40566468 # (10, 2, 1) 96
+0.40374663 # (9, 5, 0) 72
+0.40185554 # (9, 5, 1) 72
+0.39999077 # (10, 2, 2) 32
+0.39815172 # (10, 3, 0) 72
+0.39633781 # (10, 3, 1) 144
+0.39104130 # (10, 3, 2) 96
+0.38932243 # (8, 7, 1) 96
+0.38762602 # (9, 5, 3) 48
+0.38595160 # (10, 4, 0) 72
+0.38429870 # (10, 4, 1) 120
+0.38266685 # (10, 3, 3) 72
+0.37946456 # (10, 4, 2) 48
+0.37789327 # (11, 0, 0) 78
+0.37634134 # (11, 1, 0) 120
+0.37480838 # (11, 1, 1) 48
+0.37179782 # (11, 2, 0) 144
+0.37031949 # (11, 2, 1) 144
+0.36741498 # (8, 8, 0) 12
+0.36598812 # (11, 2, 2) 144
+0.36457776 # (11, 3, 0) 48
+0.36318357 # (11, 3, 1) 120
+0.36180526 # (10, 4, 4) 48
+0.36044253 # (9, 6, 4) 48
+0.35909507 # (11, 3, 2) 168
+0.35644489 # (10, 6, 0) 48
+0.35514161 # (11, 4, 0) 96
+0.35385253 # (11, 4, 1) 96
+0.35257738 # (11, 3, 3) 72
+0.35131592 # (10, 6, 2) 48
+0.35006790 # (11, 4, 2) 96
+0.34883309 # (9, 6, 5) 48
+0.34640217 # (12, 0, 0) 30
+0.34520561 # (12, 1, 0) 96
+0.34402137 # (12, 1, 1) 192
+0.34284923 # (11, 5, 1) 56
+0.34168899 # (12, 2, 0) 24
+0.34054045 # (12, 2, 1) 168
+0.33940342 # (11, 5, 2) 120
+0.33716311 # (12, 2, 2) 72
+0.33605946 # (12, 3, 0) 144
+0.33496658 # (12, 3, 1) 96
+0.33388430 # (11, 5, 3) 96
+0.33175083 # (12, 3, 2) 72
+0.33069932 # (11, 6, 1) 96
+0.32862595 # (12, 4, 0) 24
+0.32760378 # (12, 4, 1) 192
+0.32659109 # (12, 3, 3) 108
+0.32558774 # (9, 9, 1) 24
+0.32459358 # (12, 4, 2) 96
+0.32360847 # (10, 8, 1) 96
+0.32263227 # (11, 6, 3) 120
+0.32070609 # (10, 8, 2) 48
+0.31975585 # (13, 0, 0) 78
+0.31881400 # (13, 1, 0) 144
+0.31788043 # (13, 1, 1) 120
+0.31695501 # (10, 6, 6) 24
+0.31603763 # (13, 2, 0) 168
+0.31512817 # (13, 2, 1) 144
+0.31333255 # (12, 4, 4) 24
+0.31244618 # (13, 2, 2) 48
+0.31156728 # (13, 3, 0) 96
+0.31069576 # (13, 3, 1) 120
+0.30983152 # (12, 6, 0) 72
+0.30897444 # (12, 6, 1) 120
+0.30812444 # (13, 3, 2) 144
+0.30644528 # (12, 6, 2) 48
+0.30561593 # (13, 4, 0) 192
+0.30479327 # (13, 4, 1) 144
+0.30397722 # (13, 3, 3) 48
+0.30236460 # (13, 4, 2) 192
+0.30156785 # (10, 9, 3) 48
+0.29999308 # (8, 8, 8) 8
+0.29921488 # (12, 7, 0) 48
+0.29844271 # (13, 5, 0) 240
+0.29767649 # (13, 5, 1) 96
+0.29691614 # (14, 0, 0) 54
+0.29616159 # (14, 1, 0) 120
+0.29541276 # (14, 1, 1) 120
+0.29393199 # (14, 2, 0) 84
+0.29319990 # (14, 2, 1) 144
+0.29247326 # (12, 7, 3) 72
+0.29175199 # (13, 5, 3) 96
+0.29103603 # (14, 2, 2) 48
+0.29032532 # (14, 3, 0) 96
+0.28961979 # (14, 3, 1) 240
+0.28822402 # (12, 8, 0) 24
+0.28753367 # (14, 3, 2) 240
+0.28684825 # (13, 5, 4) 96
+0.28616770 # (11, 9, 3) 72
+0.28549198 # (14, 4, 0) 72
+0.28482102 # (14, 4, 1) 96
+0.28415478 # (14, 3, 3) 72
+0.28283618 # (14, 4, 2) 96
+0.28218374 # (12, 8, 3) 96
+0.28153578 # (13, 7, 0) 120
+0.28089227 # (13, 7, 1) 96
+0.27961838 # (14, 5, 0) 192
+0.27898789 # (14, 5, 1) 144
+0.27773962 # (12, 8, 4) 48
+0.27712173 # (15, 0, 0) 150
+0.27650795 # (15, 1, 0) 96
+0.27589823 # (15, 1, 1) 120
+0.27529253 # (14, 4, 4) 48
+0.27469079 # (15, 2, 0) 120
+0.27409299 # (15, 2, 1) 240
+0.27290900 # (14, 6, 0) 24
+0.27232273 # (15, 2, 2) 144
+0.27174022 # (15, 3, 0) 216
+0.27116143 # (15, 3, 1) 48
+0.27058633 # (14, 6, 2) 72
+0.27001486 # (14, 5, 4) 144
+0.26944701 # (15, 3, 2) 96
+0.26776470 # (15, 4, 0) 144
+0.26721090 # (15, 4, 1) 132
+0.26666051 # (15, 3, 3) 104
+0.26611352 # (12, 10, 0) 72
+0.26556987 # (15, 4, 2) 168
+0.26502955 # (14, 7, 1) 144
+0.26395871 # (14, 6, 4) 96
+0.26342814 # (14, 7, 2) 144
+0.26290076 # (15, 5, 0) 144
+0.26237653 # (15, 5, 1) 168
+0.26133741 # (12, 10, 3) 48
+0.26082246 # (15, 5, 2) 192
+0.25980162 # (16, 0, 0) 6
+0.25929568 # (16, 1, 0) 192
+0.25879268 # (16, 1, 1) 96
+0.25829260 # (15, 5, 3) 96
+0.25779540 # (16, 2, 0) 96
+0.25730107 # (16, 2, 1) 216
+0.25680957 # (15, 6, 1) 72
+0.25583496 # (16, 2, 2) 96
+0.25535179 # (16, 3, 0) 96
+0.25487135 # (16, 3, 1) 240
+0.25439362 # (13, 7, 7) 48
+0.25391856 # (14, 6, 6) 24
+0.25344615 # (16, 3, 2) 264
+0.25297637 # (15, 6, 3) 192
+0.25204460 # (16, 4, 0) 48
+0.25158255 # (16, 4, 1) 96
+0.25112304 # (16, 3, 3) 144
+0.25066604 # (15, 7, 1) 120
+0.25021152 # (16, 4, 2) 96
+0.24975947 # (15, 6, 4) 72
+0.24930986 # (15, 7, 2) 168
+0.24841787 # (12, 10, 6) 48
+0.24797545 # (16, 5, 0) 240
+0.24753539 # (16, 5, 1) 96
+0.24709766 # (15, 7, 3) 72
+0.24622912 # (16, 5, 2) 192
+0.24579827 # (15, 6, 5) 144
+0.24494332 # (16, 4, 4) 36
+0.24451918 # (17, 0, 0) 102
+0.24409723 # (17, 1, 0) 240
+0.24367746 # (17, 1, 1) 96
+0.24325984 # (16, 6, 0) 48
+0.24284437 # (17, 2, 0) 216
+0.24243102 # (17, 2, 1) 168
+0.24161060 # (16, 6, 2) 120
+0.24120351 # (17, 2, 2) 192
+0.24079846 # (17, 3, 0) 72
+0.24039545 # (17, 3, 1) 192
+0.23999446 # (14, 10, 2) 56
+0.23959547 # (16, 6, 3) 96
+0.23919846 # (17, 3, 2) 144
+0.23841032 # (12, 12, 4) 24
+0.23801916 # (17, 4, 0) 192
+0.23762993 # (17, 4, 1) 240
+0.23724259 # (17, 3, 3) 72
+0.23685714 # (16, 6, 4) 96
+0.23647357 # (17, 4, 2) 144
+0.23609185 # (15, 9, 2) 96
+0.23533393 # (14, 10, 4) 48
+0.23495770 # (14, 9, 6) 96
+0.23458326 # (17, 5, 0) 312
+0.23421061 # (17, 5, 1) 144
+0.23347061 # (16, 6, 5) 120
+0.23310323 # (17, 5, 2) 144
+0.23237364 # (16, 8, 0) 24
+0.23201140 # (17, 4, 4) 240
+0.23165086 # (15, 9, 4) 96
+0.23129199 # (17, 5, 3) 96
+0.23093478 # (18, 0, 0) 102
+0.23057922 # (18, 1, 0) 168
+0.23022530 # (18, 1, 1) 264
+0.22952232 # (18, 2, 0) 48
+0.22917324 # (18, 2, 1) 288
+0.22882574 # (17, 5, 4) 96
+0.22847982 # (15, 9, 5) 72
+0.22813546 # (18, 2, 2) 72
+0.22779266 # (18, 3, 0) 120
+0.22745140 # (18, 3, 1) 144
+0.22677345 # (16, 8, 4) 48
+0.22643674 # (18, 3, 2) 96
+0.22610153 # (17, 7, 0) 180
+0.22576780 # (17, 7, 1) 144
+0.22543554 # (18, 4, 0) 48
+0.22510475 # (18, 4, 1) 336
+0.22477541 # (18, 3, 3) 216
+0.22412104 # (18, 4, 2) 120
+0.22379599 # (16, 8, 5) 96
+0.22347235 # (16, 9, 3) 120
+0.22315011 # (17, 7, 3) 120
+0.22250979 # (18, 5, 0) 168
+0.22219170 # (18, 5, 1) 240
+0.22155957 # (12, 12, 8) 24
+0.22124552 # (18, 5, 2) 192
+0.22093281 # (17, 8, 1) 192
+0.22062142 # (15, 11, 3) 96
+0.22031134 # (18, 4, 4) 144
+0.22000256 # (17, 8, 2) 96
+0.21969508 # (18, 5, 3) 72
+0.21908397 # (18, 6, 0) 120
+0.21878032 # (19, 0, 0) 126
+0.21847792 # (19, 1, 0) 216
+0.21817678 # (19, 1, 1) 104
+0.21787688 # (18, 6, 2) 48
+0.21757822 # (19, 2, 0) 240
+0.21728078 # (19, 2, 1) 144
+0.21639572 # (19, 2, 2) 288
+0.21610309 # (19, 3, 0) 144
+0.21581165 # (19, 3, 1) 192
+0.21552139 # (16, 10, 4) 48
+0.21523229 # (18, 7, 0) 120
+0.21494436 # (19, 3, 2) 336
+0.21437193 # (18, 6, 4) 96
+0.21408743 # (19, 4, 0) 192
+0.21380406 # (19, 4, 1) 192
+0.21352181 # (19, 3, 3) 72
+0.21296065 # (19, 4, 2) 240
+0.21268172 # (18, 7, 3) 96
+0.21212714 # (16, 8, 8) 24
+0.21185147 # (18, 6, 5) 96
+0.21157687 # (19, 5, 0) 240
+0.21130334 # (19, 5, 1) 120
+0.21103087 # (18, 8, 0) 48
+0.21075944 # (18, 8, 1) 264
+0.21048907 # (19, 5, 2) 192
+0.20995142 # (18, 8, 2) 108
+0.20968413 # (19, 4, 4) 144
+0.20941787 # (15, 13, 0) 120
+0.20915261 # (19, 5, 3) 192
+0.20888837 # (18, 6, 6) 72
+0.20862512 # (19, 6, 0) 72
+0.20836286 # (19, 6, 1) 240
+0.20784130 # (20, 0, 0) 30
+0.20758198 # (20, 1, 0) 240
+0.20732364 # (20, 1, 1) 192
+0.20706625 # (15, 13, 3) 48
+0.20680982 # (20, 2, 0) 168
+0.20655434 # (20, 2, 1) 216
+0.20629981 # (19, 6, 3) 192
+0.20579355 # (20, 2, 2) 48
+0.20554182 # (20, 3, 0) 192
+0.20529100 # (20, 3, 1) 192
+0.20504111 # (19, 7, 1) 144
+0.20454404 # (20, 3, 2) 240
+0.20429685 # (19, 7, 2) 240
+0.20380516 # (20, 4, 0) 72
+0.20356065 # (20, 4, 1) 144
+0.20331701 # (20, 3, 3) 96
+0.20307424 # (19, 7, 3) 216
+0.20283234 # (20, 4, 2) 96
+0.20259130 # (18, 9, 4) 120
+0.20235112 # (19, 6, 5) 120
+0.20187332 # (18, 10, 0) 72
+0.20163568 # (20, 5, 0) 336
+0.20139888 # (20, 5, 1) 288
+0.20116291 # (15, 11, 9) 48
+0.20092777 # (18, 10, 2) 72
+0.20069345 # (20, 5, 2) 192
+0.20045995 # (18, 9, 5) 144
+0.19999538 # (20, 4, 4) 32
+0.19976431 # (19, 6, 6) 144
+0.19953403 # (20, 5, 3) 288
+0.19930455 # (19, 7, 5) 96
+0.19907586 # (20, 6, 0) 72
+0.19884795 # (20, 6, 1) 240
+0.19862083 # (17, 10, 7) 96
+0.19816890 # (20, 6, 2) 144
+0.19794410 # (21, 0, 0) 270
+0.19772005 # (21, 1, 0) 96
+0.19749676 # (21, 1, 1) 120
+0.19705245 # (21, 2, 0) 96
+0.19683142 # (21, 2, 1) 384
+0.19617275 # (21, 2, 2) 240
+0.19595466 # (21, 3, 0) 252
+0.19573729 # (21, 3, 1) 144
+0.19552065 # (20, 6, 4) 96
+0.19530472 # (20, 7, 2) 144
+0.19508951 # (21, 3, 2) 168
+0.19466121 # (16, 14, 2) 96
+0.19444812 # (21, 4, 0) 96
+0.19423572 # (21, 4, 1) 312
+0.19402402 # (21, 3, 3) 192
+0.19381301 # (18, 10, 6) 48
+0.19360269 # (21, 4, 2) 360
+0.19339305 # (19, 10, 1) 96
+0.19297580 # (20, 8, 0) 72
+0.19276819 # (20, 8, 1) 192
+0.19256125 # (21, 5, 0) 96
+0.19235497 # (21, 5, 1) 168
+0.19214935 # (20, 8, 2) 120
+0.19194439 # (18, 12, 1) 192
+0.19174009 # (21, 5, 2) 240
+0.19133343 # (20, 6, 6) 72
+0.19113106 # (21, 4, 4) 144
+0.19092934 # (20, 7, 5) 240
+0.19072826 # (21, 5, 3) 120
+0.19032799 # (21, 6, 0) 216
+0.19012880 # (21, 6, 1) 96
+0.18973228 # (20, 8, 4) 48
+0.18953495 # (21, 6, 2) 192
+0.18933824 # (21, 5, 4) 240
+0.18914213 # (19, 11, 1) 96
+0.18894664 # (22, 0, 0) 78
+0.18875175 # (22, 1, 0) 240
+0.18855746 # (22, 1, 1) 312
+0.18817067 # (22, 2, 0) 120
+0.18797817 # (22, 2, 1) 240
+0.18778626 # (21, 7, 0) 168
+0.18759493 # (21, 7, 1) 216
+0.18740419 # (22, 2, 2) 48
+0.18721403 # (22, 3, 0) 144
+0.18702444 # (22, 3, 1) 336
+0.18645913 # (22, 3, 2) 288
+0.18627183 # (20, 7, 7) 96
+0.18608509 # (21, 7, 3) 72
+0.18589891 # (22, 4, 0) 144
+0.18571329 # (22, 4, 1) 192
+0.18552822 # (22, 3, 3) 168
+0.18515975 # (22, 4, 2) 144
+0.18497633 # (21, 8, 0) 96
+0.18479346 # (21, 8, 1) 336
+0.18461112 # (19, 11, 5) 104
+0.18424807 # (22, 5, 0) 360
+0.18406735 # (22, 5, 1) 192
+0.18370749 # (16, 16, 0) 12
+0.18352835 # (22, 5, 2) 192
+0.18334973 # (21, 8, 3) 192
+0.18317164 # (21, 7, 5) 144
+0.18299406 # (22, 4, 4) 144
+0.18281700 # (20, 9, 6) 144
+0.18264045 # (22, 5, 3) 192
+0.18228888 # (22, 6, 0) 48
+0.18211385 # (22, 6, 1) 384
+0.18193933 # (21, 9, 0) 120
+0.18176531 # (21, 9, 1) 120
+0.18159179 # (22, 6, 2) 120
+0.18141876 # (22, 5, 4) 240
+0.18124623 # (21, 9, 2) 144
+0.18090263 # (20, 8, 8) 48
+0.18073157 # (23, 0, 0) 150
+0.18056098 # (23, 1, 0) 336
+0.18039088 # (23, 1, 1) 216
+0.18022126 # (18, 12, 8) 48
+0.18005212 # (23, 2, 0) 144
+0.17988345 # (23, 2, 1) 240
+0.17954754 # (22, 6, 4) 168
+0.17938028 # (23, 2, 2) 144
+0.17921349 # (23, 3, 0) 120
+0.17904717 # (23, 3, 1) 216
+0.17871591 # (21, 10, 0) 120
+0.17855097 # (23, 3, 2) 288
+0.17822244 # (20, 12, 0) 48
+0.17805886 # (23, 4, 0) 384
+0.17789573 # (23, 4, 1) 288
+0.17773305 # (23, 3, 3) 72
+0.17757081 # (22, 8, 0) 96
+0.17740901 # (23, 4, 2) 360
+0.17724766 # (21, 10, 3) 168
+0.17692626 # (22, 8, 2) 96
+0.17676622 # (20, 12, 3) 96
+0.17660661 # (23, 5, 0) 264
+0.17644744 # (23, 5, 1) 96
+0.17628869 # (22, 6, 6) 72
+0.17613037 # (22, 8, 3) 216
+0.17597248 # (23, 5, 2) 288
+0.17565796 # (20, 12, 4) 48
+0.17550133 # (23, 4, 4) 192
+0.17534512 # (21, 11, 0) 96
+0.17518933 # (23, 5, 3) 216
+0.17503395 # (22, 8, 4) 96
+0.17487898 # (23, 6, 0) 144
+0.17472443 # (23, 6, 1) 360
+0.17441654 # (18, 12, 10) 48
+0.17426321 # (23, 6, 2) 384
+0.17411028 # (23, 5, 4) 192
+0.17395775 # (21, 11, 3) 120
+0.17365390 # (22, 8, 5) 192
+0.17350256 # (23, 6, 3) 192
+0.17320108 # (24, 0, 0) 30
+0.17305093 # (24, 1, 0) 96
+0.17290117 # (24, 1, 1) 204
+0.17275179 # (23, 7, 1) 192
+0.17260280 # (24, 2, 0) 96
+0.17245420 # (24, 2, 1) 336
+0.17230598 # (23, 7, 2) 192
+0.17201068 # (24, 2, 2) 192
+0.17186360 # (24, 3, 0) 288
+0.17171690 # (24, 3, 1) 216
+0.17157057 # (23, 7, 3) 168
+0.17142462 # (22, 10, 2) 56
+0.17127903 # (24, 3, 2) 192
+0.17113382 # (23, 6, 5) 240
+0.17084450 # (24, 4, 0) 24
+0.17070038 # (24, 4, 1) 288
+0.17055664 # (24, 3, 3) 384
+0.17041325 # (19, 15, 3) 96
+0.17027023 # (24, 4, 2) 168
+0.17012756 # (23, 8, 2) 144
+0.16998526 # (21, 11, 6) 96
+0.16970171 # (22, 10, 4) 120
+0.16956047 # (24, 5, 0) 240
+0.16941958 # (24, 5, 1) 288
+0.16927904 # (23, 7, 5) 120
+0.16899901 # (24, 5, 2) 312
+0.16885951 # (22, 11, 1) 144
+0.16858156 # (24, 4, 4) 72
+0.16844309 # (23, 8, 4) 192
+0.16830497 # (24, 5, 3) 144
+0.16816718 # (23, 9, 1) 240
+0.16802973 # (24, 6, 0) 144
+0.16789262 # (24, 6, 1) 120
+0.16775584 # (23, 9, 2) 408
+0.16748329 # (24, 6, 2) 96
+0.16734751 # (24, 5, 4) 144
+0.16721206 # (23, 8, 5) 144
+0.16707694 # (23, 9, 3) 120
+0.16694215 # (22, 10, 6) 96
+0.16680768 # (24, 6, 3) 384
+0.16667354 # (21, 10, 9) 144
+0.16627304 # (25, 0, 0) 150
+0.16614018 # (25, 1, 0) 432
+0.16600764 # (25, 1, 1) 96
+0.16587542 # (24, 6, 4) 72
+0.16574351 # (25, 2, 0) 432
+0.16561191 # (25, 2, 1) 240
+0.16534966 # (22, 12, 2) 96
+0.16521900 # (25, 2, 2) 240
+0.16508865 # (25, 3, 0) 168
+0.16495861 # (25, 3, 1) 240
+0.16469944 # (24, 6, 5) 168
+0.16457032 # (25, 3, 2) 240
+0.16431297 # (24, 8, 0) 24
+0.16418476 # (25, 4, 0) 336
+0.16405684 # (25, 4, 1) 192
+0.16392922 # (25, 3, 3) 72
+0.16380189 # (24, 8, 2) 192
+0.16367486 # (25, 4, 2) 192
+0.16354813 # (23, 9, 6) 192
+0.16329555 # (24, 6, 6) 108
+0.16316969 # (24, 8, 3) 240
+0.16304413 # (25, 5, 0) 360
+0.16291886 # (25, 5, 1) 192
+0.16279387 # (18, 18, 2) 24
+0.16266917 # (22, 13, 0) 168
+0.16254476 # (25, 5, 2) 336
+0.16229679 # (24, 8, 4) 96
+0.16217323 # (25, 4, 4) 240
+0.16204995 # (24, 9, 1) 96
+0.16192695 # (25, 5, 3) 264
+0.16180423 # (20, 16, 2) 96
+0.16168179 # (25, 6, 0) 216
+0.16155963 # (25, 6, 1) 264
+0.16131613 # (22, 12, 6) 120
+0.16119480 # (25, 6, 2) 288
+0.16107374 # (25, 5, 4) 360
+0.16095295 # (21, 15, 1) 96
+0.16071218 # (22, 13, 4) 144
+0.16059220 # (25, 6, 3) 144
+0.16035304 # (20, 16, 4) 48
+0.16023387 # (24, 9, 4) 144
+0.16011495 # (25, 7, 0) 288
+0.15999631 # (25, 7, 1) 224
+0.15987792 # (26, 0, 0) 78
+0.15975980 # (26, 1, 0) 360
+0.15964194 # (26, 1, 1) 240
+0.15940700 # (26, 2, 0) 144
+0.15928992 # (26, 2, 1) 240
+0.15917309 # (24, 9, 5) 144
+0.15905653 # (25, 7, 3) 120
+0.15894021 # (26, 2, 2) 120
+0.15882416 # (26, 3, 0) 144
+0.15870835 # (26, 3, 1) 384
+0.15847751 # (20, 12, 12) 24
+0.15836246 # (26, 3, 2) 480
+0.15824766 # (25, 8, 1) 192
+0.15813311 # (23, 9, 9) 120
+0.15801881 # (26, 4, 0) 168
+0.15790476 # (26, 4, 1) 288
+0.15779096 # (26, 3, 3) 120
+0.15756408 # (26, 4, 2) 144
+0.15745101 # (25, 6, 6) 96
+0.15733818 # (25, 8, 3) 312
+0.15722560 # (25, 7, 5) 240
+0.15700115 # (26, 5, 0) 408
+0.15688929 # (26, 5, 1) 192
+0.15666627 # (24, 8, 8) 24
+0.15655512 # (26, 5, 2) 288
+0.15644421 # (25, 9, 0) 288
+0.15633353 # (25, 9, 1) 144
+0.15622309 # (26, 4, 4) 48
+0.15611288 # (23, 12, 6) 120
+0.15600290 # (26, 5, 3) 384
+0.15578364 # (26, 6, 0) 96
+0.15567436 # (26, 6, 1) 288
+0.15556530 # (25, 8, 5) 288
+0.15545648 # (25, 9, 3) 96
+0.15534788 # (26, 6, 2) 120
+0.15523951 # (26, 5, 4) 192
+0.15513137 # (22, 15, 3) 144
+0.15491576 # (24, 12, 0) 72
+0.15480829 # (26, 6, 3) 192
+0.15470104 # (25, 9, 4) 228
+0.15459402 # (25, 7, 7) 96
+0.15448722 # (24, 12, 2) 120
+0.15438064 # (26, 7, 0) 360
+0.15427428 # (26, 7, 1) 264
+0.15406222 # (26, 6, 4) 144
+0.15395652 # (27, 0, 0) 318
+0.15385103 # (27, 1, 0) 144
+0.15374576 # (27, 1, 1) 288
+0.15353587 # (27, 2, 0) 168
+0.15343125 # (27, 2, 1) 480
+0.15322264 # (24, 12, 4) 48
+0.15311865 # (27, 2, 2) 240
+0.15301488 # (27, 3, 0) 240
+0.15291132 # (27, 3, 1) 120
+0.15280796 # (26, 8, 0) 192
+0.15270482 # (26, 8, 1) 288
+0.15260188 # (27, 3, 2) 96
+0.15239664 # (26, 8, 2) 144
+0.15229432 # (27, 4, 0) 192
+0.15219221 # (27, 4, 1) 312
+0.15209031 # (27, 3, 3) 216
+0.15198861 # (26, 6, 6) 48
+0.15188712 # (27, 4, 2) 384
+0.15178582 # (26, 7, 5) 288
+0.15148316 # (25, 8, 8) 144
+0.15138267 # (27, 5, 0) 240
+0.15128239 # (27, 5, 1) 288
+0.15118230 # (26, 8, 4) 192
+0.15108241 # (26, 9, 0) 120
+0.15098272 # (27, 5, 2) 264
+0.15078393 # (20, 18, 6) 48
+0.15068482 # (27, 4, 4) 480
+0.15058592 # (25, 11, 4) 144
+0.15048721 # (27, 5, 3) 96
+0.15029036 # (27, 6, 0) 240
+0.15019223 # (27, 6, 1) 288
+0.14999654 # (16, 16, 16) 8
+0.14989898 # (27, 6, 2) 240
+0.14980161 # (27, 5, 4) 384
+0.14970443 # (25, 11, 5) 144
+0.14960744 # (24, 14, 0) 48
+0.14951064 # (26, 9, 4) 312
+0.14941403 # (27, 6, 3) 360
+0.14922136 # (26, 10, 0) 240
+0.14912530 # (26, 10, 1) 192
+0.14902943 # (27, 7, 0) 168
+0.14893375 # (27, 7, 1) 240
+0.14883825 # (26, 10, 2) 96
+0.14874293 # (27, 6, 4) 240
+0.14864779 # (27, 7, 2) 288
+0.14845807 # (28, 0, 0) 54
+0.14836348 # (28, 1, 0) 192
+0.14826907 # (28, 1, 1) 192
+0.14817484 # (27, 7, 3) 120
+0.14808079 # (28, 2, 0) 120
+0.14798692 # (28, 2, 1) 384
+0.14789323 # (27, 6, 5) 192
+0.14770638 # (28, 2, 2) 120
+0.14761322 # (28, 3, 0) 96
+0.14752023 # (28, 3, 1) 504
+0.14742743 # (25, 13, 1) 96
+0.14724233 # (28, 3, 2) 360
+0.14715005 # (26, 11, 1) 192
+0.14696599 # (28, 4, 0) 84
+0.14687422 # (28, 4, 1) 432
+0.14678263 # (28, 3, 3) 144
+0.14669120 # (27, 7, 5) 240
+0.14659995 # (28, 4, 2) 144
+0.14650887 # (25, 12, 6) 192
+0.14641795 # (26, 11, 3) 336
+0.14623663 # (24, 14, 6) 72
+0.14614622 # (28, 5, 0) 384
+0.14605598 # (28, 5, 1) 408
+0.14596590 # (27, 9, 1) 168
+0.14587599 # (26, 10, 6) 96
+0.14578625 # (28, 5, 2) 144
+0.14569668 # (27, 9, 2) 144
+0.14551802 # (28, 4, 4) 48
+0.14542893 # (24, 15, 4) 144
+0.14534001 # (28, 5, 3) 336
+0.14525126 # (27, 9, 3) 240
+0.14516266 # (28, 6, 0) 96
+0.14507423 # (28, 6, 1) 360
+0.14498596 # (26, 11, 5) 240
+0.14480990 # (28, 6, 2) 240
+0.14472211 # (28, 5, 4) 336
+0.14463447 # (27, 9, 4) 144
+0.14454700 # (27, 7, 7) 168
+0.14437253 # (28, 6, 3) 264
+0.14428554 # (27, 10, 1) 240
+0.14411201 # (24, 16, 0) 24
+0.14402548 # (28, 7, 0) 336
+0.14393911 # (28, 7, 1) 192
+0.14385290 # (27, 9, 5) 144
+0.14376683 # (28, 6, 4) 240
+0.14368093 # (28, 7, 2) 192
+0.14359517 # (27, 10, 3) 168
+0.14342412 # (26, 10, 8) 96
+0.14333883 # (29, 0, 0) 174
+0.14325368 # (29, 1, 0) 312
+0.14316869 # (29, 1, 1) 144
+0.14308385 # (22, 18, 6) 72
+0.14299916 # (29, 2, 0) 360
+0.14291462 # (29, 2, 1) 480
+0.14274599 # (28, 8, 0) 72
+0.14266190 # (29, 2, 2) 336
+0.14257796 # (29, 3, 0) 240
+0.14249416 # (29, 3, 1) 240
+0.14241051 # (28, 8, 2) 96
+0.14232701 # (24, 14, 9) 120
+0.14224366 # (29, 3, 2) 528
+0.14207739 # (28, 6, 6) 72
+0.14199447 # (29, 4, 0) 384
+0.14191170 # (29, 4, 1) 192
+0.14182907 # (29, 3, 3) 168
+0.14166425 # (29, 4, 2) 288
+0.14158206 # (23, 18, 3) 96
+0.14141809 # (28, 8, 4) 96
+0.14133632 # (28, 9, 0) 192
+0.14125470 # (29, 5, 0) 528
+0.14117321 # (29, 5, 1) 152
+0.14109187 # (24, 16, 6) 96
+0.14101066 # (28, 9, 2) 384
+0.14092960 # (29, 5, 2) 192
+0.14076789 # (26, 14, 0) 120
+0.14068724 # (29, 4, 4) 240
+0.14060674 # (28, 9, 3) 240
+0.14052637 # (29, 5, 3) 288
+0.14044614 # (26, 14, 2) 96
+0.14036604 # (29, 6, 0) 120
+0.14028608 # (29, 6, 1) 240
+0.14004703 # (29, 6, 2) 480
+0.13996761 # (29, 5, 4) 324
+0.13988833 # (21, 21, 1) 72
+0.13980919 # (28, 10, 0) 192
+0.13973018 # (28, 10, 1) 288
+0.13965130 # (29, 6, 3) 216
+0.13949395 # (28, 10, 2) 144
+0.13941547 # (27, 12, 4) 192
+0.13933712 # (29, 7, 0) 288
+0.13925891 # (29, 7, 1) 216
+0.13910288 # (29, 6, 4) 336
+0.13902506 # (29, 7, 2) 336
+0.13886981 # (24, 16, 8) 48
+0.13879238 # (28, 8, 7) 192
+0.13871508 # (27, 13, 0) 144
+0.13863791 # (29, 7, 3) 336
+0.13856087 # (30, 0, 0) 150
+0.13848395 # (30, 1, 0) 288
+0.13840717 # (30, 1, 1) 336
+0.13825398 # (30, 2, 0) 96
+0.13817757 # (30, 2, 1) 288
+0.13810129 # (29, 8, 1) 336
+0.13802514 # (27, 13, 3) 72
+0.13794912 # (30, 2, 2) 120
+0.13787322 # (30, 3, 0) 504
+0.13779744 # (30, 3, 1) 192
+0.13764626 # (28, 8, 8) 48
+0.13757086 # (30, 3, 2) 144
+0.13749558 # (29, 8, 3) 432
+0.13742043 # (29, 7, 5) 192
+0.13734540 # (30, 4, 0) 120
+0.13727049 # (30, 4, 1) 240
+0.13719570 # (30, 3, 3) 192
+0.13704649 # (30, 4, 2) 240
+0.13697207 # (29, 8, 4) 240
+0.13689777 # (29, 9, 0) 216
+0.13682359 # (29, 9, 1) 240
+0.13667560 # (30, 5, 0) 168
+0.13660178 # (30, 5, 1) 480
+0.13645450 # (28, 12, 0) 24
+0.13638104 # (30, 5, 2) 432
+0.13630769 # (29, 8, 5) 288
+0.13623447 # (29, 9, 3) 168
+0.13616136 # (30, 4, 4) 144
+0.13608837 # (28, 10, 7) 192
+0.13601550 # (30, 5, 3) 312
+0.13587011 # (30, 6, 0) 216
+0.13579759 # (30, 6, 1) 240
+0.13572518 # (29, 9, 4) 192
+0.13565289 # (29, 7, 7) 192
+0.13558071 # (30, 6, 2) 48
+0.13550865 # (30, 5, 4) 552
+0.13543671 # (29, 10, 1) 144
+0.13529316 # (28, 12, 4) 72
+0.13522156 # (30, 6, 3) 384
+0.13515007 # (28, 9, 9) 192
+0.13507870 # (29, 9, 5) 120
+0.13500743 # (28, 10, 8) 144
+0.13493628 # (30, 7, 0) 144
+0.13486524 # (30, 7, 1) 504
+0.13472350 # (30, 6, 4) 96
+0.13465280 # (30, 7, 2) 384
+0.13458221 # (29, 8, 7) 360
+0.13451173 # (27, 15, 1) 96
+0.13437110 # (29, 10, 4) 192
+0.13430095 # (30, 7, 3) 192
+0.13409116 # (31, 0, 0) 198
+0.13402145 # (31, 1, 0) 336
+0.13395185 # (31, 1, 1) 216
+0.13388235 # (30, 8, 0) 144
+0.13381296 # (31, 2, 0) 528
+0.13374368 # (31, 2, 1) 288
+0.13360545 # (30, 8, 2) 132
+0.13353649 # (31, 2, 2) 288
+0.13346764 # (31, 3, 0) 144
+0.13339889 # (31, 3, 1) 360
+0.13333026 # (30, 6, 6) 104
+0.13326172 # (30, 8, 3) 144
+0.13319330 # (31, 3, 2) 432
+0.13305676 # (24, 20, 0) 72
+0.13298865 # (31, 4, 0) 240
+0.13292064 # (31, 4, 1) 288
+0.13285274 # (31, 3, 3) 192
+0.13278494 # (30, 8, 4) 168
+0.13271724 # (31, 4, 2) 360
+0.13264965 # (30, 9, 1) 120
+0.13251477 # (28, 14, 2) 144
+0.13244749 # (30, 9, 2) 288
+0.13238031 # (31, 5, 0) 528
+0.13231323 # (31, 5, 1) 192
+0.13217938 # (30, 8, 5) 432
+0.13211260 # (31, 5, 2) 432
+0.13197936 # (28, 12, 8) 96
+0.13191289 # (31, 4, 4) 144
+0.13184651 # (29, 12, 3) 192
+0.13178024 # (31, 5, 3) 192
+0.13171407 # (28, 14, 4) 144
+0.13164800 # (31, 6, 0) 168
+0.13158203 # (31, 6, 1) 312
+0.13145038 # (30, 10, 0) 144
+0.13138470 # (31, 6, 2) 480
+0.13131913 # (31, 5, 4) 192
+0.13125365 # (29, 9, 9) 96
+0.13118827 # (30, 10, 2) 168
+0.13112298 # (29, 10, 8) 192
+0.13105779 # (31, 6, 3) 240
+0.13086282 # (30, 10, 3) 240
+0.13079802 # (31, 7, 0) 336
+0.13073331 # (31, 7, 1) 288
+0.13066871 # (24, 20, 6) 48
+0.13060419 # (31, 6, 4) 312
+0.13053978 # (31, 7, 2) 360
+0.13041123 # (30, 10, 4) 192
+0.13034710 # (30, 9, 6) 288
+0.13028306 # (28, 15, 3) 216
+0.13021912 # (31, 7, 3) 312
+0.13009152 # (30, 11, 0) 264
+0.13002785 # (31, 6, 5) 384
+0.12990081 # (32, 0, 0) 6
+0.12983743 # (32, 1, 0) 480
+0.12977414 # (32, 1, 1) 384
+0.12971095 # (27, 17, 3) 96
+0.12964784 # (32, 2, 0) 192
+0.12958483 # (32, 2, 1) 384
+0.12952191 # (30, 11, 3) 144
+0.12939634 # (32, 2, 2) 96
+0.12933369 # (32, 3, 0) 144
+0.12927114 # (32, 3, 1) 528
+0.12920867 # (31, 7, 5) 240
+0.12914630 # (30, 10, 6) 96
+0.12908402 # (32, 3, 2) 240
+0.12902182 # (29, 14, 1) 144
+0.12889770 # (32, 4, 0) 96
+0.12883578 # (32, 4, 1) 432
+0.12877394 # (32, 3, 3) 144
+0.12871219 # (31, 9, 1) 192
+0.12865053 # (32, 4, 2) 216
+0.12858896 # (30, 12, 1) 192
+0.12852748 # (31, 9, 2) 504
+0.12840478 # (30, 12, 2) 72
+0.12834357 # (32, 5, 0) 528
+0.12828244 # (32, 5, 1) 336
+0.12822139 # (31, 9, 3) 120
+0.12809957 # (32, 5, 2) 408
+0.12803878 # (27, 18, 1) 192
+0.12791748 # (32, 4, 4) 96
+0.12785695 # (30, 11, 6) 192
+0.12779652 # (32, 5, 3) 300
+0.12773616 # (31, 7, 7) 144
+0.12767590 # (32, 6, 0) 96
+0.12761571 # (32, 6, 1) 312
+0.12755562 # (31, 10, 1) 360
+0.12743568 # (32, 6, 2) 240
+0.12737583 # (32, 5, 4) 192
+0.12731608 # (29, 15, 0) 240
+0.12725640 # (31, 9, 5) 288
+0.12719681 # (26, 14, 14) 48
+0.12713730 # (32, 6, 3) 360
+0.12707788 # (31, 10, 3) 432
+0.12695928 # (28, 12, 12) 24
+0.12690011 # (32, 7, 0) 288
+0.12684101 # (32, 7, 1) 384
+0.12678200 # (29, 15, 3) 168
+0.12672308 # (32, 6, 4) 264
+0.12666423 # (32, 7, 2) 288
+0.12660547 # (31, 9, 6) 216
+0.12648819 # (30, 12, 6) 192
+0.12642967 # (30, 10, 9) 192
+0.12637123 # (32, 7, 3) 264
+0.12631287 # (31, 11, 1) 152
+0.12619640 # (32, 6, 5) 384
+0.12613829 # (31, 11, 2) 336
+0.12602230 # (32, 8, 0) 48
+0.12596442 # (33, 0, 0) 390
+0.12590663 # (33, 1, 0) 144
+0.12584891 # (33, 1, 1) 408
+0.12579128 # (32, 8, 2) 96
+0.12573372 # (33, 2, 0) 120
+0.12567624 # (33, 2, 1) 312
+0.12556152 # (32, 6, 6) 144
+0.12550428 # (33, 2, 2) 432
+0.12544711 # (33, 3, 0) 360
+0.12539003 # (33, 3, 1) 144
+0.12533302 # (30, 14, 2) 120
+0.12527609 # (29, 16, 2) 336
+0.12521924 # (33, 3, 2) 240
+0.12510576 # (32, 8, 4) 96
+0.12504914 # (33, 4, 0) 192
+0.12499259 # (33, 4, 1) 576
+0.12493613 # (33, 3, 3) 192
+0.12487973 # (30, 12, 8) 72
+0.12482342 # (33, 4, 2) 600
+0.12476718 # (31, 10, 7) 192
+0.12465493 # (30, 14, 4) 168
+0.12459892 # (32, 8, 5) 192
+0.12454298 # (33, 5, 0) 264
+0.12448712 # (33, 5, 1) 240
+0.12437562 # (28, 18, 3) 168
+0.12431998 # (33, 5, 2) 432
+0.12420893 # (24, 20, 12) 48
+0.12415352 # (33, 4, 4) 528
+0.12409818 # (32, 7, 7) 192
+0.12404292 # (33, 5, 3) 120
+0.12398772 # (32, 10, 0) 240
+0.12393261 # (33, 6, 0) 432
+0.12387756 # (33, 6, 1) 264
+0.12376769 # (32, 10, 2) 96
+0.12371287 # (33, 6, 2) 192
+0.12365812 # (33, 5, 4) 528
+0.12360344 # (31, 13, 1) 192
+0.12354883 # (30, 14, 6) 72
+0.12349429 # (32, 10, 3) 336
+0.12343983 # (33, 6, 3) 432
+0.12327687 # (32, 8, 7) 240
+0.12322270 # (33, 7, 0) 144
+0.12316859 # (33, 7, 1) 384
+0.12311456 # (32, 10, 4) 192
+0.12306060 # (33, 6, 4) 288
+0.12300671 # (33, 7, 2) 216
+0.12289914 # (30, 12, 10) 144
+0.12284546 # (32, 11, 0) 288
+0.12279185 # (32, 11, 1) 384
+0.12273831 # (33, 7, 3) 144
+0.12263144 # (32, 11, 2) 192
+0.12257811 # (33, 6, 5) 240
+0.12247166 # (32, 8, 8) 36
+0.12241854 # (33, 8, 0) 192
+0.12236549 # (33, 8, 1) 672
+0.12231250 # (31, 13, 5) 192
+0.12225959 # (34, 0, 0) 102
+0.12220674 # (34, 1, 0) 336
+0.12215396 # (34, 1, 1) 288
+0.12204861 # (34, 2, 0) 240
+0.12199604 # (34, 2, 1) 576
+0.12194354 # (33, 8, 3) 144
+0.12189110 # (33, 7, 5) 168
+0.12183873 # (34, 2, 2) 96
+0.12178643 # (34, 3, 0) 240
+0.12173419 # (34, 3, 1) 432
+0.12162992 # (32, 12, 0) 48
+0.12157789 # (34, 3, 2) 576
+0.12152592 # (33, 9, 0) 240
+0.12147402 # (33, 9, 1) 168
+0.12142218 # (34, 4, 0) 216
+0.12137042 # (34, 4, 1) 288
+0.12131871 # (34, 3, 3) 360
+0.12121551 # (34, 4, 2) 168
+0.12116400 # (32, 12, 3) 144
+0.12111256 # (33, 8, 5) 192
+0.12106119 # (33, 9, 3) 360
+0.12095864 # (34, 5, 0) 552
+0.12090746 # (34, 5, 1) 240
+0.12080530 # (32, 12, 4) 120
+0.12075432 # (34, 5, 2) 192
+0.12070340 # (33, 9, 4) 288
+0.12065254 # (33, 7, 7) 216
+0.12060175 # (34, 4, 4) 192
+0.12055103 # (33, 10, 0) 240
+0.12050037 # (34, 5, 3) 480
+0.12039923 # (34, 6, 0) 72
+0.12034876 # (34, 6, 1) 432
+0.12029835 # (32, 13, 1) 336
+0.12024801 # (33, 9, 5) 192
+0.12019773 # (34, 6, 2) 192
+0.12014751 # (34, 5, 4) 240
+0.12009735 # (33, 10, 3) 144
+0.11999723 # (28, 20, 4) 56
+0.11994726 # (34, 6, 3) 192
+0.11989736 # (33, 8, 7) 288
+0.11984751 # (31, 11, 11) 144
+0.11979773 # (32, 12, 6) 96
+0.11974801 # (34, 7, 0) 480
+0.11969836 # (34, 7, 1) 504
+0.11959923 # (34, 6, 4) 144
+0.11954976 # (34, 7, 2) 480
+0.11950035 # (33, 11, 0) 264
+0.11945100 # (33, 11, 1) 336
+0.11935248 # (30, 13, 12) 120
+0.11930331 # (34, 7, 3) 480
+0.11920516 # (24, 24, 8) 24
+0.11915618 # (34, 6, 5) 384
+0.11910725 # (32, 13, 5) 288
+0.11905839 # (33, 11, 3) 144
+0.11900958 # (34, 8, 0) 192
+0.11896084 # (34, 8, 1) 384
+0.11891215 # (31, 15, 6) 144
+0.11881496 # (34, 8, 2) 240
+0.11876646 # (35, 0, 0) 270
+0.11871801 # (35, 1, 0) 504
+0.11866962 # (35, 1, 1) 96
+0.11862130 # (34, 6, 6) 72
+0.11857303 # (35, 2, 0) 456
+0.11852482 # (35, 2, 1) 288
+0.11842857 # (32, 12, 8) 96
+0.11838054 # (35, 2, 2) 288
+0.11833256 # (35, 3, 0) 288
+0.11828464 # (35, 3, 1) 288
+0.11823678 # (34, 8, 4) 144
+0.11818898 # (34, 9, 0) 168
+0.11814124 # (35, 3, 2) 504
+0.11804593 # (30, 18, 4) 96
+0.11799836 # (35, 4, 0) 384
+0.11795084 # (35, 4, 1) 384
+0.11790339 # (35, 3, 3) 96
+0.11780865 # (35, 4, 2) 384
+0.11776136 # (34, 9, 3) 192
+0.11766697 # (28, 20, 8) 48
+0.11761985 # (33, 12, 4) 384
+0.11757279 # (35, 5, 0) 444
+0.11752579 # (35, 5, 1) 360
+0.11747885 # (28, 18, 12) 96
+0.11743196 # (34, 9, 4) 384
+0.11738513 # (35, 5, 2) 288
+0.11729163 # (34, 10, 0) 312
+0.11724497 # (35, 4, 4) 240
+0.11719836 # (33, 13, 0) 144
+0.11715180 # (35, 5, 3) 360
+0.11710531 # (34, 10, 2) 144
+0.11705886 # (35, 6, 0) 240
+0.11701248 # (35, 6, 1) 432
+0.11687364 # (35, 6, 2) 480
+0.11682748 # (35, 5, 4) 384
+0.11678136 # (33, 13, 3) 144
+0.11673530 # (32, 12, 10) 120
+0.11668930 # (34, 8, 7) 384
+0.11664335 # (35, 6, 3) 240
+0.11655161 # (34, 10, 4) 144
+0.11650583 # (34, 9, 6) 240
+0.11646009 # (35, 7, 0) 504
+0.11641441 # (35, 7, 1) 240
+0.11632322 # (35, 6, 4) 408
+0.11627770 # (35, 7, 2) 240
+0.11618682 # (32, 16, 0) 24
+0.11614146 # (34, 11, 2) 288
+0.11609615 # (33, 12, 7) 144
+0.11605090 # (35, 7, 3) 264
+0.11600570 # (34, 8, 8) 240
+0.11596055 # (33, 14, 0) 144
+0.11591546 # (35, 6, 5) 696
+0.11582543 # (30, 18, 8) 96
+0.11578049 # (35, 8, 0) 432
+0.11573561 # (35, 8, 1) 192
+0.11569077 # (33, 11, 9) 216
+0.11564599 # (34, 10, 6) 96
+0.11560126 # (35, 8, 2) 288
+0.11555659 # (33, 14, 3) 336
+0.11546739 # (36, 0, 0) 102
+0.11542287 # (36, 1, 0) 144
+0.11537840 # (36, 1, 1) 384
+0.11533398 # (35, 7, 5) 192
+0.11528961 # (36, 2, 0) 168
+0.11524529 # (36, 2, 1) 600
+0.11520103 # (34, 11, 5) 192
+0.11511265 # (36, 2, 2) 264
+0.11506854 # (36, 3, 0) 480
+0.11502447 # (36, 3, 1) 216
+0.11498046 # (35, 9, 1) 264
+0.11489259 # (36, 3, 2) 288
+0.11484873 # (35, 9, 2) 432
+0.11476116 # (36, 4, 0) 48
+0.11471745 # (36, 4, 1) 288
+0.11467379 # (36, 3, 3) 576
+0.11463018 # (35, 9, 3) 144
+0.11458662 # (36, 4, 2) 288
+0.11454311 # (32, 17, 2) 240
+0.11449965 # (34, 9, 9) 120
+0.11441287 # (34, 10, 8) 96
+0.11436956 # (36, 5, 0) 288
+0.11432629 # (36, 5, 1) 504
+0.11428308 # (35, 7, 7) 224
+0.11423991 # (30, 18, 10) 72
+0.11419679 # (36, 5, 2) 504
+0.11415372 # (35, 10, 1) 480
+0.11406773 # (36, 4, 4) 72
+0.11402481 # (35, 10, 2) 432
+0.11398193 # (36, 5, 3) 288
+0.11393911 # (35, 9, 5) 288
+0.11389633 # (36, 6, 0) 120
+0.11385360 # (36, 6, 1) 240
+0.11381092 # (35, 10, 3) 384
+0.11372570 # (36, 6, 2) 144
+0.11368316 # (36, 5, 4) 288
+0.11364067 # (35, 8, 7) 240
+0.11359823 # (33, 15, 5) 192
+0.11351348 # (36, 6, 3) 504
+0.11347118 # (35, 9, 6) 240
+0.11338672 # (32, 16, 8) 48
+0.11334457 # (36, 7, 0) 192
+0.11330245 # (36, 7, 1) 336
+0.11326039 # (35, 11, 1) 144
+0.11321837 # (36, 6, 4) 96
+0.11317640 # (36, 7, 2) 672
+0.11313447 # (35, 11, 2) 480
+0.11305076 # (34, 14, 0) 180
+0.11300898 # (35, 8, 8) 192
+0.11296724 # (36, 7, 3) 264
+0.11292555 # (35, 11, 3) 288
+0.11288390 # (34, 14, 2) 144
+0.11284230 # (36, 6, 5) 192
+0.11280074 # (34, 11, 9) 288
+0.11271777 # (36, 8, 0) 48
+0.11267635 # (36, 8, 1) 720
+0.11263498 # (35, 11, 4) 288
+0.11259366 # (33, 15, 7) 144
+0.11255237 # (36, 8, 2) 336
+0.11251114 # (31, 20, 2) 192
+0.11246995 # (33, 14, 9) 216
+0.11238770 # (36, 6, 6) 216
+0.11234665 # (37, 0, 0) 222
+0.11230564 # (37, 1, 0) 528
+0.11226467 # (37, 1, 1) 288
+0.11218288 # (37, 2, 0) 216
+0.11214205 # (37, 2, 1) 336
+0.11206052 # (36, 8, 4) 120
+0.11201982 # (37, 2, 2) 432
+0.11197917 # (37, 3, 0) 240
+0.11193856 # (37, 3, 1) 384
+0.11189799 # (32, 16, 10) 96
+0.11185747 # (36, 9, 2) 312
+0.11181700 # (37, 3, 2) 456
+0.11173618 # (32, 18, 6) 120
+0.11169583 # (37, 4, 0) 576
+0.11165553 # (37, 4, 1) 480
+0.11161527 # (37, 3, 3) 96
+0.11157506 # (34, 14, 6) 120
+0.11153488 # (37, 4, 2) 336
+0.11149476 # (34, 15, 3) 240
+0.11137463 # (36, 9, 4) 192
+0.11133468 # (37, 5, 0) 576
+0.11129477 # (37, 5, 1) 288
+0.11125490 # (36, 10, 0) 168
+0.11121507 # (36, 10, 1) 288
+0.11117529 # (37, 5, 2) 240
+0.11109585 # (36, 10, 2) 240
+0.11105619 # (37, 4, 4) 336
+0.11101658 # (36, 9, 5) 168
+0.11097701 # (37, 5, 3) 336
+0.11089799 # (37, 6, 0) 288
+0.11085855 # (37, 6, 1) 528
+0.11077979 # (24, 24, 16) 24
+0.11074047 # (37, 6, 2) 432
+0.11070119 # (37, 5, 4) 384
+0.11066196 # (31, 21, 3) 96
+0.11062276 # (36, 10, 4) 192
+0.11058361 # (36, 9, 6) 360
+0.11054450 # (37, 6, 3) 336
+0.11046640 # (34, 16, 2) 192
+0.11042742 # (36, 11, 0) 192
+0.11038847 # (37, 7, 0) 408
+0.11034957 # (37, 7, 1) 288
+0.11031071 # (30, 22, 6) 96
+0.11027189 # (37, 6, 4) 648
+0.11023311 # (37, 7, 2) 288
+0.11015567 # (36, 8, 8) 144
+0.11011701 # (35, 14, 2) 336
+0.11007839 # (36, 11, 3) 384
+0.11003982 # (37, 7, 3) 360
+0.11000128 # (34, 16, 4) 96
+0.10996279 # (33, 18, 4) 264
+0.10992433 # (37, 6, 5) 384
+0.10984754 # (36, 10, 6) 72
+0.10980921 # (37, 8, 0) 432
+0.10977091 # (37, 8, 1) 384
+0.10973266 # (33, 15, 11) 96
+0.10965627 # (37, 8, 2) 288
+0.10961813 # (33, 18, 5) 192
+0.10954198 # (36, 12, 0) 120
+0.10950397 # (37, 6, 6) 336
+0.10946599 # (37, 8, 3) 288
+0.10942806 # (37, 7, 5) 192
+0.10939016 # (38, 0, 0) 126
+0.10935230 # (38, 1, 0) 456
+0.10931448 # (38, 1, 1) 384
+0.10923896 # (38, 2, 0) 216
+0.10920126 # (38, 2, 1) 576
+0.10916360 # (37, 9, 0) 168
+0.10912598 # (37, 9, 1) 312
+0.10908839 # (38, 2, 2) 104
+0.10905085 # (38, 3, 0) 168
+0.10901334 # (38, 3, 1) 720
+0.10893844 # (36, 12, 4) 48
+0.10890105 # (38, 3, 2) 288
+0.10886370 # (37, 8, 5) 324
+0.10882638 # (37, 9, 3) 264
+0.10878911 # (38, 4, 0) 240
+0.10875187 # (38, 4, 1) 336
+0.10871467 # (38, 3, 3) 192
+0.10864039 # (38, 4, 2) 144
+0.10860330 # (36, 13, 0) 192
+0.10856626 # (37, 9, 4) 696
+0.10852925 # (37, 7, 7) 120
+0.10845534 # (38, 5, 0) 672
+0.10841845 # (38, 5, 1) 432
+0.10830798 # (38, 5, 2) 336
+0.10827124 # (36, 13, 3) 192
+0.10823453 # (37, 9, 5) 360
+0.10819786 # (38, 4, 4) 288
+0.10816123 # (36, 10, 9) 192
+0.10812463 # (38, 5, 3) 360
+0.10805155 # (38, 6, 0) 144
+0.10801506 # (38, 6, 1) 624
+0.10797861 # (37, 8, 7) 288
+0.10794220 # (33, 15, 13) 168
+0.10790583 # (38, 6, 2) 192
+0.10786949 # (38, 5, 4) 384
+0.10783319 # (37, 9, 6) 240
+0.10776069 # (32, 20, 8) 48
+0.10772450 # (38, 6, 3) 240
+0.10768835 # (37, 11, 0) 432
+0.10765223 # (37, 11, 1) 288
+0.10761615 # (36, 14, 0) 120
+0.10758010 # (38, 7, 0) 264
+0.10754409 # (38, 7, 1) 600
+0.10747218 # (38, 6, 4) 336
+0.10743628 # (38, 7, 2) 240
+0.10740041 # (36, 11, 9) 192
+0.10736458 # (37, 11, 3) 312
+0.10729303 # (36, 14, 3) 288
+0.10725730 # (38, 7, 3) 288
+0.10718597 # (36, 12, 8) 96
+0.10715035 # (38, 6, 5) 480
+0.10711477 # (37, 11, 4) 288
+0.10707923 # (29, 21, 15) 96
+0.10704372 # (38, 8, 0) 192
+0.10700824 # (38, 8, 1) 480
+0.10697280 # (33, 15, 14) 192
+0.10690203 # (38, 8, 2) 192
+0.10686670 # (37, 12, 0) 192
+0.10683140 # (37, 12, 1) 600
+0.10679613 # (37, 11, 5) 288
+0.10676090 # (38, 6, 6) 72
+0.10672571 # (38, 8, 3) 576
+0.10669055 # (38, 7, 5) 288
+0.10658528 # (39, 0, 0) 390
+0.10655026 # (39, 1, 0) 240
+0.10651528 # (39, 1, 1) 168
+0.10648032 # (38, 8, 4) 240
+0.10644541 # (39, 2, 0) 360
+0.10641052 # (39, 2, 1) 576
+0.10634086 # (36, 14, 6) 96
+0.10630608 # (39, 2, 2) 624
+0.10627133 # (39, 3, 0) 432
+0.10623662 # (39, 3, 1) 264
+0.10616730 # (38, 8, 5) 288
+0.10613269 # (39, 3, 2) 240
+0.10606357 # (32, 16, 16) 24
+0.10602906 # (39, 4, 0) 192
+0.10599458 # (39, 4, 1) 528
+0.10596014 # (39, 3, 3) 408
+0.10592573 # (36, 12, 10) 96
+0.10589136 # (39, 4, 2) 432
+0.10585702 # (38, 7, 7) 240
+0.10578844 # (38, 10, 0) 240
+0.10575420 # (38, 10, 1) 288
+0.10571999 # (39, 5, 0) 408
+0.10568581 # (39, 5, 1) 288
+0.10565167 # (38, 10, 2) 120
+0.10561756 # (37, 12, 6) 216
+0.10558349 # (39, 5, 2) 672
+0.10551543 # (36, 16, 0) 48
+0.10548146 # (39, 4, 4) 480
+0.10544751 # (37, 13, 4) 288
+0.10541360 # (39, 5, 3) 96
+0.10537972 # (36, 16, 2) 264
+0.10534588 # (39, 6, 0) 504
+0.10531206 # (39, 6, 1) 144
+0.10524453 # (38, 10, 4) 192
+0.10521082 # (39, 6, 2) 384
+0.10517713 # (39, 5, 4) 336
+0.10514348 # (37, 13, 5) 144
+0.10507628 # (38, 11, 0) 336
+0.10504272 # (39, 6, 3) 576
+0.10497571 # (36, 16, 4) 108
+0.10494225 # (38, 11, 2) 336
+0.10490882 # (39, 7, 0) 240
+0.10487543 # (39, 7, 1) 408
+0.10484207 # (38, 8, 8) 144
+0.10480874 # (39, 6, 4) 264
+0.10477544 # (39, 7, 2) 648
+0.10470893 # (30, 26, 0) 120
+0.10467573 # (37, 12, 8) 336
+0.10464256 # (35, 17, 8) 192
+0.10460942 # (39, 7, 3) 216
+0.10457631 # (38, 10, 6) 192
+0.10454323 # (38, 11, 4) 480
+0.10451018 # (39, 6, 5) 192
+0.10444418 # (36, 12, 12) 72
+0.10441123 # (39, 8, 0) 288
+0.10437831 # (39, 8, 1) 528
+0.10434542 # (37, 13, 7) 200
+0.10431256 # (38, 12, 0) 72
+0.10427973 # (39, 8, 2) 624
+0.10424693 # (38, 11, 5) 384
+0.10418143 # (38, 12, 2) 240
+0.10414873 # (39, 6, 6) 192
+0.10411605 # (39, 8, 3) 408
+0.10408341 # (39, 7, 5) 384
+0.10401821 # (38, 12, 3) 168
+0.10398566 # (37, 15, 2) 384
+0.10392065 # (40, 0, 0) 30
+0.10388819 # (40, 1, 0) 672
+0.10385576 # (40, 1, 1) 480
+0.10382336 # (39, 9, 1) 144
+0.10379099 # (40, 2, 0) 240
+0.10375865 # (40, 2, 1) 192
+0.10372635 # (39, 9, 2) 336
+0.10366182 # (40, 2, 2) 192
+0.10362960 # (40, 3, 0) 336
+0.10359741 # (40, 3, 1) 384
+0.10356525 # (39, 9, 3) 360
+0.10353313 # (30, 26, 6) 48
+0.10350103 # (40, 3, 2) 504
+0.10346896 # (38, 13, 1) 336
+0.10340491 # (40, 4, 0) 168
+0.10337293 # (40, 4, 1) 336
+0.10334098 # (40, 3, 3) 144
+0.10330906 # (39, 7, 7) 360
+0.10327717 # (40, 4, 2) 216
+0.10324531 # (39, 10, 0) 216
+0.10321348 # (39, 10, 1) 360
+0.10314991 # (38, 12, 6) 192
+0.10311816 # (40, 5, 0) 576
+0.10308645 # (40, 5, 1) 240
+0.10305476 # (39, 9, 5) 168
+0.10299148 # (40, 5, 2) 600
+0.10295988 # (39, 10, 3) 336
+0.10289678 # (40, 4, 4) 48
+0.10286527 # (36, 16, 9) 192
+0.10283378 # (40, 5, 3) 768
+0.10280233 # (35, 19, 7) 192
+0.10277091 # (40, 6, 0) 192
+0.10273951 # (40, 6, 1) 456
+0.10270815 # (39, 9, 6) 432
+0.10264550 # (40, 6, 2) 192
+0.10261422 # (40, 5, 4) 528
+0.10258297 # (39, 11, 0) 168
+0.10255175 # (39, 11, 1) 240
+0.10252055 # (38, 14, 2) 144
+0.10248939 # (40, 6, 3) 192
+0.10245825 # (39, 11, 2) 528
+0.10236501 # (40, 7, 0) 576
+0.10233398 # (40, 7, 1) 480
+0.10230299 # (39, 11, 3) 192
+0.10227202 # (40, 6, 4) 240
+0.10224108 # (40, 7, 2) 192
+0.10221017 # (33, 23, 6) 264
+0.10214843 # (38, 14, 4) 240
+0.10211760 # (39, 10, 6) 192
+0.10208680 # (40, 7, 3) 504
+0.10205603 # (37, 17, 1) 192
+0.10199456 # (40, 6, 5) 576
+0.10196388 # (38, 13, 7) 240
+0.10190258 # (40, 8, 0) 72
+0.10187198 # (40, 8, 1) 576
+0.10184140 # (39, 12, 1) 336
+0.10181085 # (39, 11, 5) 312
+0.10178032 # (40, 8, 2) 144
+0.10174983 # (39, 12, 2) 312
+0.10171936 # (39, 10, 7) 336
+0.10165850 # (40, 6, 6) 96
+0.10162812 # (40, 8, 3) 384
+0.10159776 # (40, 7, 5) 576
+0.10156742 # (37, 15, 9) 168
+0.10153712 # (38, 14, 6) 216
+0.10150684 # (38, 13, 8) 192
+0.10147659 # (39, 11, 6) 240
+0.10141617 # (40, 8, 4) 96
+0.10138600 # (41, 0, 0) 246
+0.10135586 # (41, 1, 0) 372
+0.10132574 # (41, 1, 1) 240
+0.10129565 # (36, 18, 8) 120
+0.10126559 # (41, 2, 0) 624
+0.10123555 # (41, 2, 1) 528
+0.10117556 # (38, 12, 10) 120
+0.10114561 # (41, 2, 2) 432
+0.10111568 # (41, 3, 0) 312
+0.10108577 # (41, 3, 1) 432
+0.10102605 # (37, 18, 0) 264
+0.10099623 # (41, 3, 2) 624
+0.10093666 # (36, 20, 0) 72
+0.10090691 # (41, 4, 0) 336
+0.10087720 # (41, 4, 1) 192
+0.10084750 # (41, 3, 3) 264
+0.10081784 # (40, 10, 0) 336
+0.10078820 # (41, 4, 2) 624
+0.10075859 # (39, 10, 9) 240
+0.10069944 # (40, 10, 2) 288
+0.10066990 # (38, 15, 6) 192
+0.10064039 # (41, 5, 0) 696
+0.10061091 # (41, 5, 1) 240
+0.10058145 # (30, 22, 18) 48
+0.10055202 # (40, 10, 3) 504
+0.10052262 # (41, 5, 2) 240
+0.10046388 # (36, 20, 4) 72
+0.10043456 # (41, 4, 4) 432
+0.10040525 # (39, 12, 7) 240
+0.10037598 # (41, 5, 3) 384
+0.10034673 # (40, 10, 4) 192
+0.10031750 # (41, 6, 0) 192
+0.10028830 # (41, 6, 1) 552
+0.10022998 # (36, 18, 10) 144
+0.10020085 # (41, 6, 2) 624
+0.10017175 # (41, 5, 4) 288
+0.10014268 # (39, 11, 9) 120
+0.10008461 # (40, 11, 2) 480
+0.10005561 # (41, 6, 3) 288
diff --git a/src/geomtools/powder/calibrants/__init__.py b/src/geomtools/powder/calibrants/__init__.py
new file mode 100644
index 0000000..cbbc627
--- /dev/null
+++ b/src/geomtools/powder/calibrants/__init__.py
@@ -0,0 +1,16 @@
+from importlib.resources import files
+
+import numpy as np
+
+calibrants = files("geomtools.powder.calibrants")
+
+
+def get_calibrant(name, inverse=True):
+ fn = calibrants.joinpath(f"{name}.D")
+ d = np.loadtxt(str(fn)) / 10 # Ã… -> nm
+ return 1. / d if inverse else d
+
+
+def calibrant_list():
+ return [f.stem for f in calibrants.iterdir()
+ if f.suffix == ".D"]
diff --git a/src/geomtools/powder/find_center.py b/src/geomtools/powder/find_center.py
new file mode 100644
index 0000000..1df6428
--- /dev/null
+++ b/src/geomtools/powder/find_center.py
@@ -0,0 +1,124 @@
+import os
+from argparse import ArgumentParser
+
+import h5py
+import matplotlib.pyplot as plt
+from extra.components import AGIPD1MQuadrantMotors
+from extra_data import open_run
+from extra_geom import AGIPD_1MGeometry, JUNGFRAUGeometry
+from extra_geom.motors import AGIPD_1MMotors
+from . import PowderDiffraction
+from .calibrants import calibrant_list, get_calibrant
+from .misc import agipd_borders, jungfrau_borders
+
+DETECTOR_TYPE = {
+ "agipd1m": AGIPD_1MGeometry,
+ "jf4m": JUNGFRAUGeometry,
+}
+
+BORDER_MASK = {
+ "agipd1m": agipd_borders,
+ "jf4m": jungfrau_borders,
+}
+
+CALIBRANTS = calibrant_list()
+
+
+def guess_detector_type(detector_id):
+ for detector_type in DETECTOR_TYPE.keys():
+ if detector_type.upper() in detector_id:
+ return detector_type
+
+ raise ValueError("Unknown detector type")
+
+
+def main(argv=None):
+ parser = ArgumentParser(
+ description="The program finds the beam position by "
+ "powder diffraction data")
+
+ parser.add_argument("-g", "--geometry", required=True,
+ help="Crystfel geometry file")
+ parser.add_argument("-l", "--layout", choices=DETECTOR_TYPE.keys(),
+ help="Detector layout")
+ parser.add_argument("-c", "--calibrant", required=True,
+ choices=CALIBRANTS, help="Calibrant")
+ parser.add_argument("-o", "--output", default=".",
+ help="Output directory")
+ parser.add_argument("powder_data",
+ help="Averaged powder diffraction data")
+
+ args = parser.parse_args(argv)
+
+ with h5py.File(args.powder_data, "r") as f:
+ # conditions
+ clen = f["powderSum/conditions/cameraLen"][()]
+ photon_en = f["powderSum/conditions/photonEnergy"][()]
+ lmd = f["powderSum/conditions/waveLength"][()]
+ detector_id = f["powderSum/conditions/detectorId"][()].decode()
+ runno = f["powderSum/conditions/run"][()]
+ propno = f["powderSum/conditions/proposal"][()]
+
+ # average image
+ img_mean = f["powderSum/image/mean"][:]
+ num_frames = f["powderSum/image/numFrames"][:]
+ mask = f["powderSum/image/mask"][:].astype(bool)
+ modules = f["powderSum/image/modules"][:]
+
+ num_modules = len(modules)
+
+ print("Proposal:", propno)
+ print("Run:", runno)
+
+ if args.layout is None:
+ detector_type = guess_detector_type(detector_id)
+ else:
+ detector_type = args.layout
+
+ print("Detector Id:", detector_id)
+ print("Detector type:", detector_type)
+
+ print("Photon energy", photon_en, "(keV)")
+ print("Wave length", lmd, "(nm)")
+ print("Camera length", clen, "(m)")
+
+ geom_class = DETECTOR_TYPE[detector_type]
+ ref_geom = geom_class.from_crystfel_geom(args.geometry)
+
+ if detector_type == "agipd1m":
+ run = open_run(propno, runno, data="raw")
+ tracker = AGIPD_1MMotors(ref_geom)
+ motors = AGIPD1MQuadrantMotors(run)
+ geom = tracker.geom_at(motors.most_frequent_positions())
+ else:
+ geom = ref_geom
+
+ rings = get_calibrant(args.calibrant)
+ borders = BORDER_MASK[detector_type](num_modules)
+ pw = PowderDiffraction(geom, rings, clen, lmd,
+ border_mask=~borders, make_shadow_mask=True)
+ pw.fit(img_mean, mask, num_frames)
+ pw.refinement_info()
+
+ geom_new = pw.transform_geometry()
+ refined_geom_fn = os.path.join(
+ args.output, f"powder_center_p{propno:06d}_r{runno:04d}.geom")
+ geom_new.write_crystfel_geom(
+ refined_geom_fn,
+ clen=geom_new.metadata["clen"],
+ photon_energy=geom_new.metadata["photon_energy"],
+ adu_per_ev=geom_new.metadata["photon_energy"] / 1000,
+ )
+
+ fig, (ax1, ax2) = plt.subplots(
+ 2, 1, figsize=(12, 16), tight_layout=True,
+ gridspec_kw={'height_ratios': [3, 1]})
+ pw.draw(ax=ax1)
+ pw.draw_radial_profile(img_mean, mask, ax=ax2)
+ report_img_fn = os.path.join(
+ args.output, f"powder_center_p{propno:06d}_r{runno:04d}.png")
+ fig.savefig(report_img_fn, dpi=300)
+
+
+if __name__ == "__main__":
+ main()
diff --git a/src/geomtools/powder/image_tools.py b/src/geomtools/powder/image_tools.py
new file mode 100644
index 0000000..a492c87
--- /dev/null
+++ b/src/geomtools/powder/image_tools.py
@@ -0,0 +1,100 @@
+import numpy as np
+from scipy.ndimage import gaussian_filter
+from skimage import measure
+
+from .misc import azint_fast, pixel_radii
+
+
+def make_shadow_mask(img, bad_pixels=None, threshold=0.5,
+ max_noise_spot_size=100):
+ """Makes shadow mask."""
+ max_bad_spot_size = 10000
+ if bad_pixels is None:
+ bad_pixels = np.zeros_like(img, dtype=bool)
+ mask = img > threshold
+ for modno in range(len(mask)):
+ mod_labels = measure.label(mask[modno], background=0)
+ area = np.bincount(mod_labels.ravel())
+
+ for ix in range(1, len(area)):
+ if area[ix] < max_noise_spot_size:
+ mask[modno][mod_labels == ix] = False
+
+ mod_labels = measure.label(mask[modno], background=1)
+ area = np.bincount(mod_labels.ravel())
+
+ for ix in range(1, len(area)):
+ if area[ix] < max_noise_spot_size:
+ mask[modno][mod_labels == ix] = True
+ elif area[ix] < max_bad_spot_size:
+ flag = mod_labels == ix
+ mask[modno][flag] = ~bad_pixels[modno][flag]
+
+ return ~mask
+
+
+def threshold_gaussian(geom, image, bad_pixels=None, shadow=None,
+ sigma=7, snr=1, return_background=False):
+ """Binarize image."""
+ if bad_pixels is None:
+ bad_pixels = np.zeros_like(image, dtype=bool)
+ if shadow is None:
+ shadow = np.zeros_like(image, dtype=bool)
+
+ ri = pixel_radii(geom, dtype=int)
+ v = azint_fast(ri, image, shadow | bad_pixels)
+
+ mask = bad_pixels & ~shadow
+
+ img = np.copy(image)
+ img[mask] = v[ri[mask]]
+ img[shadow] = 0
+
+ shadow_ext = (
+ gaussian_filter(shadow.astype(float), sigma, axes=(1, 2)) > 0.01)
+
+ lg_im = np.log(img)
+ minv = np.min(lg_im[~shadow_ext])
+ lg_im[~np.isfinite(lg_im)] = minv
+
+ mu = gaussian_filter(lg_im, sigma, axes=(1, 2))
+ sig = gaussian_filter(lg_im * lg_im, sigma, axes=(1, 2))
+ sig = np.sqrt(sig - mu * mu)
+
+ fg = lg_im > mu + snr * sig
+ fg = fg & ~(bad_pixels | shadow_ext)
+
+ return (fg, mu, sig) if return_background else fg
+
+
+def find_rings(r, img, peaks, min_area=30):
+ r_min = peaks[0] / 2
+
+ labels = np.zeros_like(img, dtype=int)
+ num_images = len(img)
+ for modno in range(num_images):
+ mod_labels = labels[modno]
+ mod_labels[:] = measure.label(img[modno], background=0)
+ area = np.bincount(mod_labels.ravel())
+
+ for ix in range(1, len(area)):
+ if area[ix] < min_area:
+ mod_labels[mod_labels == ix] = False
+
+ for label in np.unique(mod_labels):
+ if label == 0:
+ continue
+
+ px = mod_labels == label
+ rr = r[modno][px]
+
+ r_mean = np.mean(rr)
+ dr = np.abs(peaks - r_mean)
+
+ peak_ix = np.argmin(dr)
+ if dr[peak_ix] > r_min:
+ mod_labels[px] = 0
+ else:
+ mod_labels[px] = peak_ix + 1
+
+ return labels
diff --git a/src/geomtools/powder/misc.py b/src/geomtools/powder/misc.py
new file mode 100644
index 0000000..76a7e9e
--- /dev/null
+++ b/src/geomtools/powder/misc.py
@@ -0,0 +1,69 @@
+import numpy as np
+
+
+def azint_fast(ri, img, mask=None, nan=np.nan):
+ if mask is None:
+ r, v = ri, img
+ else:
+ r = ri[~mask]
+ v = img[~mask]
+ n = np.bincount(r)
+ rdf = np.bincount(r, weights=v)
+ np.divide(rdf, n, where=n > 0, out=rdf)
+ rdf[n == 0] = nan
+ return rdf
+
+
+def pixel_radii(geom, dtype=None):
+ xyz = geom.get_pixel_positions() / geom.pixel_size
+ x = xyz[..., 0]
+ y = xyz[..., 1]
+ r = np.sqrt(x * x + y * y)
+ if dtype is not None:
+ r = r.astype(dtype)
+ return r
+
+
+def reciprocal_to_dectector(s, lmd, clen):
+ # s = q / (2pi)
+ # theta = 2. * np.arcsin(lmd * s / 2.) # theta = 2arcsin(lmd * q / 4pi)
+ # r = (clen / pixel_size) * np.tan(theta) # r = L * tan(theta)
+ lmd_s = lmd * s
+ lmd_s_2 = lmd_s * lmd_s
+ r = clen * lmd_s * np.sqrt(4.0 - lmd_s_2) / (2.0 - lmd_s_2)
+ return r
+
+
+def agipd_borders(num_modules=16):
+ borders = np.zeros([num_modules, 512, 128], bool)
+ borders[:, 0::64, :] = True
+ borders[:, 63::64, :] = True
+ borders[:, :, 0] = True
+ borders[:, :, 127] = True
+ return borders
+
+
+def agipd_pixel_area():
+ px_area = np.ones([512, 128], np.float32)
+ px_area[63::64, :] = 2.0
+ px_area[64::64, :] = 2.0
+ return px_area
+
+
+def jungfrau_borders(num_modules=8):
+ borders = np.zeros([num_modules, 512, 1024], bool)
+ borders[:, :, 0::256] = True
+ borders[:, :, 255::256] = True
+ borders[:, 0::256, :] = True
+ borders[:, 255::256, :] = True
+ return borders
+
+
+def jungfrau_pixel_area():
+ fs = np.ones(512, np.float32)
+ fs[255::256] = 2.0
+ fs[256::256] = 2.0
+ ss = np.ones(1024, np.float32)
+ ss[255::256] = 2.0
+ ss[256::256] = 2.0
+ return np.outer(fs, ss)
diff --git a/src/geomtools/powder/powder.py b/src/geomtools/powder/powder.py
new file mode 100644
index 0000000..af7cb67
--- /dev/null
+++ b/src/geomtools/powder/powder.py
@@ -0,0 +1,146 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from pyFAI.azimuthalIntegrator import AzimuthalIntegrator
+
+from .image_tools import find_rings, make_shadow_mask, threshold_gaussian
+from .misc import pixel_radii, reciprocal_to_dectector
+from .refine import fit_rings, refine_clen
+
+
+class PowderDiffraction:
+ def __init__(self, geom, peaks, clen, wave_length,
+ min_area=30, sigma=7, snr=1, border_mask=None,
+ make_shadow_mask=True):
+ self.geom = geom
+ self.peaks = peaks
+ self.clen = clen
+ self.lmd = wave_length
+ self.min_area = min_area
+ self.sigma = sigma
+ self.snr = snr
+ self.make_shadow_mask = make_shadow_mask
+ self.labels = None
+ self.beam_center = None
+ self.shadow = None
+ if border_mask is None:
+ border_mask = np.ones(geom.expected_data_shape, bool)
+ self.border_mask = border_mask
+
+ def fit(self, img, mask, num_images=1):
+ if self.make_shadow_mask:
+ threshold = 1.0 / np.sqrt(num_images)
+ if np.size(threshold) > 1:
+ slice_thr = (slice(None),) + (None,) * (img.ndim - 1)
+ threshold = threshold[slice_thr]
+ self.shadow = make_shadow_mask(img, ~mask, threshold)
+
+ fg = threshold_gaussian(
+ self.geom, img, mask, self.shadow, self.sigma, self.snr)
+
+ fg = fg & self.border_mask
+
+ xyz = self.geom.get_pixel_positions()
+ ri = np.linalg.norm(xyz[..., :2], axis=-1)
+
+ peaks = reciprocal_to_dectector(self.peaks, self.lmd, self.clen)
+ self.labels = find_rings(ri, fg, peaks, self.min_area)
+ res = fit_rings(xyz[..., :2], self.labels, peaks)
+
+ self.beam_center = (float(res[0]), float(res[1]))
+ self.peak_no = res[2]
+ self.peak_radii = res[3]
+
+ self.clen_hat = refine_clen(
+ self.peaks, self.lmd, xyz[..., :2], self.beam_center, self.labels
+ )
+ self.lmd_hat = self.lmd
+ self.peak_radii /= self.geom.pixel_size
+
+ def transform_geometry(self):
+ if self.beam_center is None:
+ raise ValueError("Rings are not fitted yet.")
+
+ xc, yc = self.beam_center
+ geom = self.geom.offset([-xc, -yc])
+
+ hc_e = 1.2398419843320026
+ geom.metadata.update(
+ {
+ "clen": self.clen_hat,
+ "photon_energy": hc_e / self.lmd_hat * 1000,
+ }
+ )
+ motor_positions = getattr(self, "motor_positions", None)
+ if motor_positions is not None:
+ geom.motor_positions = motor_positions
+ return geom
+
+ def draw(self, ax=None):
+ geom = self.transform_geometry()
+
+ cm = plt.cm.Set1.copy()
+ cm.set_under([0.3, 0.3, 0.3])
+ cm.set_over([0, 0, 0])
+ lab_draw = np.copy(self.labels)
+ bg = lab_draw == 0
+ lab_draw = (lab_draw - 1) % cm.N + 1
+ lab_draw[bg] = cm.N + 1
+ ax = geom.plot_data(lab_draw, ax=ax, colorbar=False, cmap=cm,
+ interpolation="none", vmin=1, vmax=cm.N)
+
+ peaks = reciprocal_to_dectector(
+ self.peaks[self.peak_no - 1], self.lmd_hat,
+ self.clen_hat / geom.pixel_size)
+
+ for r in peaks:
+ ax.add_patch(
+ plt.Circle((0, 0), r, fill=False, ls=":", lw=1, ec="w"))
+
+ return ax
+
+ def refinement_info(self):
+ if self.beam_center is None:
+ raise ValueError("Rings are not fitted yet.")
+ hc_e = 1.2398419843320026
+ xc, yc = np.array(self.beam_center) / self.geom.pixel_size
+ print(f"Beam center: ({xc:.2f}, {yc:.2f}) px")
+ print(f"Camera len: {self.clen:.4f} -> {self.clen_hat:.4f} (m)")
+ print(f"Photon energy: {hc_e / self.lmd_hat:.3f} eV")
+
+ def draw_radial_profile(self, img, mask, ax=None, figsize=(10, 5)):
+ geom = self.transform_geometry()
+
+ peaks = reciprocal_to_dectector(
+ self.peaks, self.lmd_hat, self.clen_hat * 1000
+ ) # in mm
+ peaks = peaks[:max(self.peak_no)]
+
+ ai = AzimuthalIntegrator(
+ detector=geom.to_pyfai_detector(),
+ dist=self.clen_hat, # sample-detector distance (m)
+ wavelength=self.lmd_hat * 1e-10, # wavelength (m)
+ )
+
+ rmin = np.min(pixel_radii(geom)[~mask]) * geom.pixel_size * 1000
+
+ msk = mask
+ if self.shadow is not None:
+ msk |= self.shadow
+
+ rint, intens = ai.integrate1d(
+ img.reshape(-1, 128),
+ mask=msk.reshape(-1, 128),
+ npt=1000,
+ unit="r_mm",
+ radial_range=(rmin, peaks[-1] + 1),
+ )
+ if ax is None:
+ fig, ax = plt.subplots(1, 1, figsize=figsize)
+
+ for r in peaks:
+ ax.axvline(r, c="gray", lw=0.5, ls=":")
+
+ ax.semilogy(rint, intens)
+ ax.set_xlim(rmin, peaks[-1] + 1)
+
+ return ax
diff --git a/src/geomtools/powder/powdersum.py b/src/geomtools/powder/powdersum.py
new file mode 100644
index 0000000..bea57e6
--- /dev/null
+++ b/src/geomtools/powder/powdersum.py
@@ -0,0 +1,356 @@
+import multiprocessing as mp
+import os
+import re
+import time
+from argparse import ArgumentParser
+
+import h5py
+import numpy as np
+import psutil
+from extra_data import RunDirectory, open_run
+from extra_data.read_machinery import find_proposal
+
+from . import shmemarray as shmem
+from .misc import agipd_pixel_area, jungfrau_pixel_area
+
+SASE = {
+ "SA1": "SA1", "FXE": "SA1", "SPB": "SA1",
+ "SA2": "SA1", "HED": "SA2", "MID": "SA2",
+ "SA3": "SA3", "SCS": "SA3", "SQS": "SA3", "SXP": "SA3",
+}
+
+XGM = {
+ "SA1": "SA1_XTD2_XGM/XGM/DOOCS",
+ "SA2": "",
+ "SA3": "",
+}
+
+ZMOTOR = {
+ "SPB_DET_AGIPD1M-1": (
+ "SPB_IRU_AGIPD1M1/MDL/DATA_SELECTOR",
+ "spbIruAgipd1MMotorZStepper.actualPosition",
+ ),
+}
+MIN_CLEN = {"SPB_DET_AGIPD1M-1": 0.118}
+
+PX_AREA = {
+ "agipd1m": agipd_pixel_area,
+ "jf4m": jungfrau_pixel_area,
+}
+
+
+class ImageAgg:
+ def __init__(
+ self,
+ nproc,
+ det_source,
+ shape,
+ adu_per_unit=1,
+ rounding_threshold=None,
+ px_area=None,
+ mask=None,
+ ):
+ self.det_source = det_source
+ self.adu_per_unit = adu_per_unit
+ if rounding_threshold is None:
+ self.rounding = False
+ self.rounding_shift = 0.0
+ else:
+ self.rounding = True
+ self.rounding_shift = rounding_threshold - 0.5
+ self.shape = shape
+ self.nproc = nproc
+ if mask is None:
+ mask = np.zeros(shape, int)
+
+ self.mask = mask == 0
+ arr_shape = (nproc,) + tuple(shape)
+ if px_area is None:
+ px_area = np.ones(shape, np.float32)
+
+ self.px_area = px_area
+
+ self.part_mean = shmem.empty(arr_shape, float)
+ self.part_deviation = shmem.empty(arr_shape, float)
+ self.part_count = shmem.empty(arr_shape, int)
+ self.part_num_frames = shmem.empty(nproc, int)
+
+ def compute(self, args):
+ i, dc = args
+ mean = self.part_mean[i]
+ mean[:] = 0
+ deviation = self.part_deviation[i]
+ deviation[:] = 0
+ count = self.part_count[i]
+ count[:] = 0
+ self.part_num_frames[i] = np.sum(
+ dc[self.det_source].data_counts(index_group="image",
+ labelled=False)
+ )
+
+ for tid, data in dc.trains():
+ a = data[self.det_source]["image.data"] / self.adu_per_unit
+ msk = (data[self.det_source]["image.mask"] == 0) & self.mask
+
+ a[~msk] = 0.0
+ if self.rounding:
+ msk &= a > -self.rounding_shift - 0.5
+ np.round(a - self.rounding_shift, out=a)
+ np.clip(a, 0.0, None, out=a)
+
+ mean += np.sum(a, 0)
+ deviation += np.sum(a * a, 0)
+ count += np.sum(msk.astype(int), 0)
+ return i
+
+ def finish(self):
+ self.r = type("Moments", (), {})
+ self.count = np.sum(self.part_count, 0)
+ nz = self.count > 0
+
+ self.mean = np.zeros(self.shape, np.float32)
+ self.mean = np.divide(
+ np.sum(self.part_mean, 0), self.count,
+ out=self.mean, where=nz
+ )
+ self.mean /= self.px_area
+
+ self.deviation = np.zeros(self.shape, np.float32)
+ self.deviation = np.divide(
+ np.sum(self.part_deviation, 0), self.count,
+ out=self.deviation, where=nz
+ )
+ self.deviation = np.sqrt(self.deviation - self.mean * self.mean)
+ self.deviation /= self.px_area
+
+ self.num_frames = np.sum(self.part_num_frames)
+
+
+def main(argv=None):
+ parser = ArgumentParser(
+ description="The program sum up images over a given run")
+ parser.add_argument("-n", "--num-proc", type=int,
+ help="The number of processes")
+ parser.add_argument("-p", "--proposal", type=int, required=True,
+ help="proposal no.")
+ parser.add_argument("-r", "--run", type=int, required=True, help="run no.")
+ parser.add_argument("-D", "--images-dir",
+ help="Directory with calibrated detector images "
+ "if different")
+ parser.add_argument("--legacy", action="store_true",
+ help="Use legacy naming for corrected detector "
+ "sources")
+ parser.add_argument("-X", "--xgm-id", help="XGM device Id")
+ parser.add_argument("--z-motor-id", help="Detector z-axis motor Id")
+ parser.add_argument("--z-motor-key", default="activePosition",
+ help="Detector z-axis motor Id")
+ parser.add_argument("-s", "--spacing-len", type=float, default=0.0,
+ help="Length of additional spacing tubes (m)")
+ parser.add_argument("-L", "--camera-len", type=float,
+ help="Sample-to-detector distance (m)")
+ parser.add_argument("-R", "--round-to-photons", action="store_true",
+ help="Convert and round intensities to photons")
+ parser.add_argument("-t", "--rounding-threshold", type=float, default=0.7,
+ help="The rounding threshold")
+ parser.add_argument("-P", "--photon-energy", type=float,
+ help="The rounding threshold")
+ parser.add_argument("-M", "--modules",
+ type=lambda s: eval(f"np.r_[{s}]").tolist(),
+ help="List of modules")
+ parser.add_argument("-l", "--layout", choices=PX_AREA.keys(),
+ help="Detector layout")
+ parser.add_argument("-m", "--mask", help="Bad pixel mask filename")
+ parser.add_argument("--mask-path", default="/entry_1/data_1/mask",
+ help="Data source path in h5 file",)
+ parser.add_argument("-o", "--output-dir",
+ help="Directory to write computed images")
+ parser.add_argument("detector_id", help="Detector Id")
+
+ args = parser.parse_args(argv)
+
+ propdir = find_proposal(f"p{args.proposal:06d}")
+
+ if args.num_proc is None:
+ num_proc = min(psutil.cpu_count(logical=False), 32)
+
+ detector_id = args.detector_id
+ inst, _, _ = detector_id.partition("_")
+ print("Instrument:", inst)
+
+ run = open_run(args.proposal, args.run, data="all")
+
+ # photon energy and wave length
+ hc_e = 1.2398419843320026
+ photon_en = args.photon_energy
+ if photon_en:
+ # wave length
+ lmd = hc_e / photon_en
+ else:
+ xgm_id = args.xgm_id
+ if not xgm_id:
+ xgm_id = XGM[SASE[inst]]
+ print("XGM:", xgm_id)
+
+ # wave length
+ lmd = run[xgm_id, "pulseEnergy.wavelengthUsed"].as_single_value()
+ # photon energy
+ photon_en = hc_e / lmd
+
+ # camera length
+ clen = args.camera_len
+ if clen is None:
+ # read from motors
+ z_motor_id = args.z_motor_id
+ if z_motor_id:
+ z_motor_key = args.z_motor_key
+ else:
+ z_motor_id, z_motor_key = ZMOTOR[detector_id]
+
+ print(f"Detector Z motor: {z_motor_id}.{z_motor_key}")
+
+ z_motor = run[z_motor_id]
+ clen0 = MIN_CLEN.get(detector_id)
+ clen = 1e-3 * z_motor.run_value(z_motor_key) + clen0 + args.spacing_len
+
+ print()
+ print(f"Photon energy: {photon_en:g} (keV)")
+ print(f"Wave length: {lmd:.4g} (nm)")
+ print(f"Camera length: {clen:.3f} (m)")
+
+ if args.images_dir:
+ images_dir = args.images_dir
+ if images_dir[0] == "^":
+ images_dir = os.path.join(propdir, images_dir[1:])
+ try_add_rundir = os.path.join(images_dir, f"r{args.run:04d}")
+ if os.path.isdir(try_add_rundir):
+ images_dir = try_add_rundir
+ dc_img = RunDirectory(images_dir)
+ else:
+ dc_img = run
+
+ # detector sources
+ source_pattern = "{detector_id}/CORR/{modno}CH0:output"
+ key_pattern = "image.*"
+
+ det_source_pattern = source_pattern.format(
+ detector_id=detector_id, modno="*")
+ if args.legacy:
+ det_source_pattern = det_source_pattern.replace("/CORR/", "/DET/")
+ det_source_pattern = det_source_pattern.replace(":output", ":xtdf")
+
+ dc_img = dc_img.select([(det_source_pattern, key_pattern)])
+
+ modules = args.modules
+ if modules is None:
+ re_pattern = re.compile(
+ "^" +
+ source_pattern.format(
+ detector_id=detector_id, modno=r"(?P<modno>\d+)") +
+ "$"
+ )
+ matches = (re_pattern.match(s) for s in dc_img.all_sources)
+ modules = sorted(int(m["modno"]) for m in matches if m is not None)
+ num_modules = len(modules)
+
+ if args.mask:
+ mask_fn = args.mask
+ if mask_fn[0] == "^":
+ mask_fn = os.path.join(propdir, mask_fn[1:])
+
+ with h5py.File(mask_fn, "r") as f:
+ mask0 = f[args.mask_path][:]
+ else:
+ mask_fn = None
+ mask0 = [None] * num_modules
+
+ print()
+ print("Detector source pattern:", det_source_pattern)
+ print("Num modules:", num_modules)
+ print("Modules:", modules)
+ print("Mask:", mask_fn)
+ print()
+ print("Num proc:", num_proc)
+
+ if args.round_to_photons:
+ adu_per_unit = photon_en
+ rounding_threshold = args.rounding_threshold
+ print("Rounding threshold:", rounding_threshold)
+ else:
+ adu_per_unit = 1
+ rounding_threshold = None
+
+ pw_aggs = []
+ for modi, modno in enumerate(modules):
+ tm0 = time.perf_counter()
+ print(f"{modno: 3d}: ", end="")
+
+ source_id = source_pattern.format(detector_id=detector_id, modno=modno)
+ dc = dc_img.select(
+ [(source_id, "image.data"), (source_id, "image.mask")],
+ require_all=True
+ )
+
+ shape = dc[source_id, "image.data"].shape[1:]
+ if args.layout is None:
+ px_area = None
+ else:
+ px_area = PX_AREA[args.layout]()
+
+ agg = ImageAgg(
+ num_proc,
+ source_id,
+ shape,
+ adu_per_unit,
+ rounding_threshold,
+ px_area=px_area,
+ mask=mask0[modi],
+ )
+
+ with mp.Pool(num_proc) as pool:
+ result_iter = pool.imap_unordered(
+ agg.compute, enumerate(dc.split_trains(num_proc))
+ )
+
+ for r in result_iter:
+ pass
+
+ agg.finish()
+ pw_aggs.append(agg)
+
+ tm1 = time.perf_counter()
+ print(f"{tm1 - tm0:.1f} s")
+
+ output_fn = f"powder_sum_p{args.proposal:06d}_r{args.run:04d}.h5"
+ if args.output_dir:
+ output_fn = os.path.join(args.output_dir, output_fn)
+ with h5py.File(output_fn, "w") as f:
+ # conditions
+ f["powderSum/conditions/cameraLen"] = clen
+ f["powderSum/conditions/photonEnergy"] = photon_en
+ f["powderSum/conditions/waveLength"] = lmd
+ f["powderSum/conditions/detectorId"] = detector_id.encode("ascii")
+ f["powderSum/conditions/run"] = args.run
+ f["powderSum/conditions/proposal"] = args.proposal
+
+ # average image
+ f["powderSum/image/mean"] = np.stack(
+ [pw_aggs[i].mean for i in range(num_modules)]
+ )
+ f["powderSum/image/std"] = np.stack(
+ [pw_aggs[i].deviation for i in range(num_modules)]
+ )
+ f["powderSum/image/count"] = np.stack(
+ [pw_aggs[i].count for i in range(num_modules)]
+ )
+ f["powderSum/image/numFrames"] = np.stack(
+ [pw_aggs[i].num_frames for i in range(num_modules)]
+ )
+ f["powderSum/image/mask"] = np.stack(
+ [~pw_aggs[i].mask | (pw_aggs[i].count == 0)
+ for i in range(num_modules)]
+ ).astype(np.uint8)
+ f["powderSum/image/modules"] = modules
+
+
+if __name__ == "__main__":
+ main()
diff --git a/src/geomtools/powder/refine.py b/src/geomtools/powder/refine.py
new file mode 100644
index 0000000..e171fe3
--- /dev/null
+++ b/src/geomtools/powder/refine.py
@@ -0,0 +1,75 @@
+import numpy as np
+
+from .misc import reciprocal_to_dectector
+
+
+def q_moments(x, y):
+ n = x.size
+
+ xx = x * x
+ yy = y * y
+ xy = x * y
+ f = xx + yy
+
+ Mx = np.sum(x)
+ My = np.sum(y)
+ Mf = np.sum(f)
+ Mxx = np.sum(xx)
+ Mxy = np.sum(xy)
+ Myy = np.sum(yy)
+ Mff = np.sum(f * f)
+ Mfx = np.sum(f * x)
+ Mfy = np.sum(f * y)
+
+ return type("PointMoments", (), dict(
+ n=n, x=Mx, y=My, f=Mf, xx=Mxx, yy=Myy,
+ ff=Mff, xy=Mxy, fx=Mfx, fy=Mfy
+ ))
+
+
+def fit_rings(xy, labels, peaks):
+ x = xy[..., 0]
+ y = xy[..., 1]
+
+ found_peaks = np.unique(labels)
+ n = len(found_peaks) + 1
+ W = np.zeros([n, n], dtype=float)
+ E = np.zeros(n, dtype=float)
+
+ for i in range(n - 2):
+ peak_no = found_peaks[i + 1]
+
+ pixels = labels == peak_no
+ M = q_moments(x[pixels], y[pixels])
+
+ j = 2 + i
+ W[i, :2] = [2.0 * M.x, 2.0 * M.y]
+ W[-2:, j] = [M.x, M.y]
+ W[i, j] = M.n
+ W[-2:, :2] += [[2 * M.xx, 2 * M.xy], [2 * M.xy, 2 * M.yy]]
+
+ E[i] = M.f
+ E[-2:] += [M.fx, M.fy]
+
+ # WU = E
+ U = np.linalg.solve(W, E)
+ xc, yc = U[0], U[1]
+
+ peak_r_hat = np.sqrt(U[2:] + xc * xc + yc * yc)
+ peak_no = found_peaks[1:]
+
+ return xc, yc, peak_no, peak_r_hat
+
+
+def refine_clen(s, lmd, xy, center, labels):
+ xc, yc = center
+ x = xy[..., 0] - xc
+ y = xy[..., 1] - yc
+
+ m = labels > 0
+ r_hat_2 = x[m] * x[m] + y[m] * y[m]
+ ix = labels[m] - 1
+
+ r = reciprocal_to_dectector(s[ix], lmd, 1.0)
+ clen_hat = np.sqrt(np.sum(r_hat_2) / np.sum(r * r))
+ return clen_hat
diff --git a/src/geomtools/powder/shmemarray.py b/src/geomtools/powder/shmemarray.py
new file mode 100644
index 0000000..a03e603
--- /dev/null
+++ b/src/geomtools/powder/shmemarray.py
@@ -0,0 +1,144 @@
+import mmap
+
+import numpy
+from packaging.version import Version
+
+NUMPY_LT_2_0 = not (Version(numpy.__version__) < Version("2.0.dev"))
+
+
+def empty_like(array, dtype=None):
+ """Create a shared memory array from the shape of array."""
+ array = numpy.asarray(array)
+ if dtype is None:
+ dtype = array.dtype
+ return anonymousmemmap(array.shape, dtype)
+
+
+def empty(shape, dtype="f8"):
+ """Create an empty shared memory array."""
+ return anonymousmemmap(shape, dtype)
+
+
+def full_like(array, value, dtype=None):
+ """Create a shared memory array with the same shape and type
+ as a given array, filled with `value`."""
+ shared = empty_like(array, dtype)
+ shared[:] = value
+ return shared
+
+
+def full(shape, value, dtype="f8"):
+ """Create a shared memory array of given shape and type,
+ filled with `value`."""
+ shared = empty(shape, dtype)
+ shared[:] = value
+ return shared
+
+
+def copy(a):
+ """Copy an array to the shared memory.
+
+ Notes
+ -----
+ copy is not always necessary because the private memory is
+ always copy-on-write.
+
+ Use :code:`a = copy(a)` to immediately dereference the old 'a'
+ on private memory
+ """
+ shared = anonymousmemmap(a.shape, dtype=a.dtype)
+ shared[:] = a[:]
+ return shared
+
+
+def fromiter(iter, dtype, count=None):
+ return copy(numpy.fromiter(iter, dtype, count))
+
+
+try:
+ # numpy >= 1.16
+ _unpickle_ctypes_type = numpy.ctypeslib.as_ctypes_type(numpy.dtype("|u1"))
+except Exception:
+ # older version numpy < 1.16
+ _unpickle_ctypes_type = numpy.ctypeslib._typecodes["|u1"]
+
+
+def __unpickle__(ai, dtype):
+ dtype = numpy.dtype(dtype)
+ tp = _unpickle_ctypes_type
+
+ # The following lb ub logic replaces ctypeslib.as_ctypes()
+ # that does not support strides.
+ if ai["strides"]:
+ lb = 0
+ ub = dtype.itemsize
+ for s, t in zip(ai["strides"], ai["shape"]):
+ if t == 0: # there is no data if any shape is 0.
+ lb, ub = 0, 0
+ break
+ if s < 0:
+ lb = lb + (t - 1) * s
+ else:
+ ub = ub + (t - 1) * s
+ else:
+ lb = 0
+ ub = dtype.itemsize
+ for s in ai["shape"]:
+ ub = ub * s
+
+ # grab a flat char array at the sharemem address,
+ # covering the memory region ai refers to.
+ ra = (tp * (ub - lb)).from_address(ai["data"][0] + lb)
+
+ # view it as what it should look like. here we assume numpy
+ # will not do crazy things like modifying gaps due to striding.
+ shm = numpy.ndarray(
+ buffer=ra, dtype=dtype, offset=-lb,
+ strides=ai["strides"], shape=ai["shape"]
+ ).view(type=anonymousmemmap)
+ return shm
+
+
+class anonymousmemmap(numpy.memmap):
+ """Arrays allocated on shared memory.
+
+ The array is stored in an anonymous memory map
+ that is shared between child-processes.
+
+ """
+
+ def __new__(subtype, shape, dtype=numpy.uint8, order="C"):
+
+ descr = numpy.dtype(dtype)
+ _dbytes = descr.itemsize
+
+ shape = numpy.atleast_1d(shape)
+ size = 1
+ for k in shape:
+ size *= k
+
+ bytes = int(size * _dbytes)
+
+ if bytes > 0:
+ mm = mmap.mmap(-1, bytes)
+ else:
+ mm = numpy.empty(0, dtype=descr)
+ self = numpy.ndarray.__new__(
+ subtype, shape, dtype=descr, buffer=mm, order=order
+ )
+ self._mmap = mm
+ return self
+
+ def __array_wrap__(self, outarr, context=None, return_scalar=False):
+ # after ufunc this won't be on shm!
+ if NUMPY_LT_2_0:
+ return numpy.ndarray.__array_wrap__(
+ self.view(numpy.ndarray), outarr, context
+ )
+ else:
+ return numpy.ndarray.__array_wrap__(
+ self.view(numpy.ndarray), outarr, context, return_scalar
+ )
+
+ def __reduce__(self):
+ return __unpickle__, (self.__array_interface__, self.dtype)
--
GitLab