diff --git a/src/geomtools/__init__.py b/src/geomtools/__init__.py
index b0abc1d6b185fb1e00920a0113acfcb4b8b6bc41..cc94950e98db73327f6e9091f1d56a7d32bc4b92 100644
--- a/src/geomtools/__init__.py
+++ b/src/geomtools/__init__.py
@@ -1,6 +1,6 @@
 # flake8: noqa E401
 from .sfx import (
     parse_crystfel_streamfile, read_crystfel_streamfile, extract_geometry,
-    plot_center_shift, plot_cell_parameters, plot_peakogram, spacing,
-    gauss2d_fit, ellipse, parse_xwiz_summary,
+    plot_center_shift, plot_cell_parameters, plot_peakogram, plot_powder,
+    spacing, gauss2d_fit, ellipse, parse_xwiz_summary, get_peak_position,
 )
diff --git a/src/geomtools/sfx/__init__.py b/src/geomtools/sfx/__init__.py
index fcfa1c0e453c4370b0461efbd76356f1a2cf4f57..76fa95bff471e297ede008cdfd0160367c240133 100644
--- a/src/geomtools/sfx/__init__.py
+++ b/src/geomtools/sfx/__init__.py
@@ -1,8 +1,10 @@
 # flake8: noqa E401
 from .crystfelio import (
-    parse_crystfel_streamfile, read_crystfel_streamfile, extract_geometry
+    parse_crystfel_streamfile, read_crystfel_streamfile, extract_geometry,
+)
+from .draw import (
+    plot_center_shift, plot_cell_parameters, plot_peakogram, plot_powder,
 )
-from .draw import plot_center_shift, plot_cell_parameters, plot_peakogram
 from .lattice import spacing
-from .misc import gauss2d_fit, ellipse
+from .misc import gauss2d_fit, ellipse, get_peak_position
 from .xwizio import parse_xwiz_summary
diff --git a/src/geomtools/sfx/crystfelio.py b/src/geomtools/sfx/crystfelio.py
index c2f0e5a5fb2e15664a17485f9f8aa9588052d225..02d530ecfce040d66b274b5af203d442d575a318 100644
--- a/src/geomtools/sfx/crystfelio.py
+++ b/src/geomtools/sfx/crystfelio.py
@@ -90,6 +90,7 @@ def parse_crystfel_streamfile(stream_filename, begin=0, end=None):
     peaks = pd.concat(peaks, ignore_index=True)
     reflexes = pd.concat(reflexes, ignore_index=True)
     lattices = pd.DataFrame(lattices)
+
     return peaks, lattices, reflexes, frame_ix
 
 
diff --git a/src/geomtools/sfx/draw.py b/src/geomtools/sfx/draw.py
index a354e21d918574b0be5a9f2fc4007277bb4b956f..c9bd27794293929f2a1e2a9caf0d3a3289abb977 100644
--- a/src/geomtools/sfx/draw.py
+++ b/src/geomtools/sfx/draw.py
@@ -99,3 +99,30 @@ def plot_peakogram(reflex, figsize=(12, 11), frameon=False, **kwargs):
     cax.set_ylabel('counts', fontsize=12)
 
     return fig, ax
+
+
+def plot_powder(peaks, counts=True, figwidth=14, frameon=False, **kwargs):
+    x = peaks.x
+    y = peaks.y
+    xbin = int(np.round(x.max() - x.min()))
+    ybin = int(np.round(y.max() - y.min()))
+    H, yedges, xedges = np.histogram2d(
+        y.values, x.values, bins=(ybin, xbin),
+        weights=None if counts else peaks.intensity
+    )
+    low, med, high = np.percentile(H.ravel(), [2, 50, 98])
+
+    height, width = H.shape
+    figsize = (figwidth, figwidth / width * height)
+    fig, ax = plt.subplots(1, 1, figsize=figsize, frameon=frameon, **kwargs)
+
+    extent = (xedges[-1], xedges[0], yedges[0], yedges[-1])
+    ax.imshow(np.flip(H, axis=1), origin='lower',
+              extent=extent, vmin=max(0, low), vmax=high)
+    ax.plot([-50, 50], [0, 0], 'r')
+    ax.plot([0, 0], [-50, 50], 'r')
+
+    ax.set_aspect('equal')
+    ax.axis(False)
+
+    return fig, ax
diff --git a/src/geomtools/sfx/misc.py b/src/geomtools/sfx/misc.py
index f5fa79e60542c7e7bf4349f80fed857d81ae3ae3..6d76cbf962a2384a286ea013ddb1c81e86e5286a 100644
--- a/src/geomtools/sfx/misc.py
+++ b/src/geomtools/sfx/misc.py
@@ -1,4 +1,5 @@
 import numpy as np
+import pandas as pd
 
 
 def gauss2d_fit(x, y):
@@ -37,3 +38,16 @@ def ellipse(p):
     x, y = u / d, Ds / d
 
     return np.array([b, a]), np.array([[x, y], [y, -x]])
+
+
+def get_peak_position(peaks, panels):
+    panels_slice = panels[['cnx', 'cny', 'fsx', 'fsy', 'ssx', 'ssy',
+                           'orig_min_fs', 'orig_min_ss']]
+    pp = peaks[['Panel', 'fs/px', 'ss/px']].join(panels_slice, on='Panel')
+    x = ((pp['fs/px'] - pp.orig_min_fs) * pp.fsx +
+         (pp['ss/px'] - pp.orig_min_ss) * pp.ssx +
+         pp.cnx)
+    y = ((pp['fs/px'] - pp.orig_min_fs) * pp.fsy +
+         (pp['ss/px'] - pp.orig_min_ss) * pp.ssy +
+         pp.cny)
+    return pd.DataFrame({'x': x, 'y': y, 'intensity': peaks.Intensity})
diff --git a/src/geomtools/sfx/report.ipynb b/src/geomtools/sfx/report.ipynb
index 5e2f5501db7569f2fb441d41155f5078e7053f84..7894afe55e4732f5868cec3740ab1cf1dfc5a78e 100644
--- a/src/geomtools/sfx/report.ipynb
+++ b/src/geomtools/sfx/report.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "id": "3852d9ec-1143-4846-9f48-11751729d547",
    "metadata": {
     "tags": [
@@ -15,7 +15,8 @@
     "stream_file = \"\"\n",
     "summary_file = \"\"\n",
     "prefix=prefix = \"\"\n",
-    "output_dir = \"\""
+    "output_dir = \"\"\n",
+    "geometry_file = \"\""
    ]
   },
   {
@@ -27,9 +28,12 @@
    "source": [
     "import os\n",
     "import matplotlib.pyplot as plt\n",
+    "import pandas as pd\n",
+    "from cfelpyutils.geometry import load_crystfel_geometry\n",
     "from geomtools import (\n",
     "    read_crystfel_streamfile, extract_geometry, parse_xwiz_summary,\n",
-    "    plot_center_shift, plot_cell_parameters, plot_peakogram, \n",
+    "    plot_center_shift, plot_cell_parameters, plot_peakogram,\n",
+    "    plot_powder, get_peak_position,\n",
     ")"
    ]
   },
@@ -44,6 +48,16 @@
     "propno = int(xwiz_conf['data']['proposal'])\n",
     "runs = eval(xwiz_conf['data']['runs'])\n",
     "\n",
+    "geom = load_crystfel_geometry(geometry_file)\n",
+    "panels = pd.DataFrame(geom.detector['panels'].values(),\n",
+    "                      index=geom.detector['panels'].keys())\n",
+    "panels.index.name='Panel'\n",
+    "\n",
+    "photon_energy = geom.beam['photon_energy']\n",
+    "clen = panels.clen.unique()\n",
+    "if len(clen) == 1:\n",
+    "    clen = clen[0]\n",
+    "\n",
     "pk, la, re, nfrm = read_crystfel_streamfile(stream_file, disp=False)"
    ]
   },
@@ -65,6 +79,9 @@
     "print(f\"Proposal: {propno}\")\n",
     "print(f\"Runs: {runs}\")\n",
     "print()\n",
+    "print(f\"Photon energy: {photon_energy} eV\")\n",
+    "print(f\"Distance to sample: {clen} m\")\n",
+    "print()\n",
     "print(summary)\n",
     "print()\n",
     "print(\"From Crystfel stream file\")\n",
@@ -127,6 +144,25 @@
     "fig, ax = plot_peakogram(re)\n",
     "plt.show()"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7f05317e",
+   "metadata": {},
+   "source": [
+    "# Powder sum"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3cdf4100",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "peak_pos = get_peak_position(pk, panels)\n",
+    "fig, ax = plot_powder(peak_pos, counts=True)"
+   ]
   }
  ],
  "metadata": {
diff --git a/src/geomtools/sfx/report.py b/src/geomtools/sfx/report.py
index af20d6a1e329a53f6b19d8efcb01ecd21419d03e..14e35bbf2887a24639b556cf862be7b3edf7039c 100644
--- a/src/geomtools/sfx/report.py
+++ b/src/geomtools/sfx/report.py
@@ -117,11 +117,7 @@ def push_geometry():
             print("Target directory exists, use another 'tag'. Quit...")
             exit(-1)
 
-        geom = extract_geometry(stream_file)
-        if geom is None:
-            print("Geometry is not found in crystfel stream file. "
-                  "Use file from xwix.summary.")
-            geom = open(xwiz_conf['geom']['file_path'], 'r')
+        geom_file = output_dir / (geom_id + ".geom")
 
         dc = open_run(propno, runs[0], include="*-DA??-*")
         for runno in runs[1:]:
@@ -138,11 +134,6 @@ def push_geometry():
             'date': date,
             'motors': motor_pos,
         }
-
-        geom_file = output_dir / (geom_id + ".geom")
-        geom_file.write_text(geom.read())
-        geom.close()
-
         geom_file.with_suffix(".yaml").write_text(yaml.dump(meta))
 
         summary_copy = (output_dir / f"{geom_id}_xwiz.summary")
@@ -150,17 +141,32 @@ def push_geometry():
     else:
         print("Dry run, only generate report...")
         output_dir = pathlib.Path.cwd()
+        geom_file = output_dir / (geom_id + ".geom")
+
+    geom = extract_geometry(stream_file)
+    if geom is None:
+        print("Geometry is not found in crystfel stream file. "
+              "Use file from xwix.summary.")
+        geom = open(xwiz_conf['geom']['file_path'], 'r')
+
+    geom_file.write_text(geom.read())
+    geom.close()
 
     notebook_filename = output_dir / f"{geom_id}_report.ipynb"
     report_parameters = dict(
         xwiz_dir=str(xwiz_dir.absolute()),
         stream_file=str(stream_file.absolute()),
         summary_file=str(summary_file.absolute()),
+        geometry_file=str(geom_file.absolute()),
         prefix=prefix,
         output_dir=str(output_dir.absolute()),
     )
     make_report(notebook_filename, report_parameters)
 
+    if args.report_only:
+        print("Clean up files...")
+        geom_file.unlink()
+
 
 if __name__ == "__main__":
     push_geometry()