diff --git a/.gitignore b/.gitignore index 6725208..1b2514f 100644 --- a/.gitignore +++ b/.gitignore @@ -248,6 +248,14 @@ test/k4FWCoreTest/**/*.root test/inputFiles/*.slcio test/gaudi_opts/testConverterConstants.py -# Files produced during running examples -*root -*png +# Files produced during running examples. Anchored on the dot so that names that +# merely end in "root" or "png" (and directories called root/) stay tracked. +*.root +*.png + +# Small deterministic input that keeps the EDM4hep notebook runnable after a fresh clone. +!DD4hepTutorials/data/ +!DD4hepTutorials/data/simplecalo2_sample.root + +# geoWebDisplay writes this next to the geometry it is showing. +viewer.cxx diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 6706137..94d84c7 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -36,4 +36,4 @@ repos: "-o", "Key4hep-Project", "-n", "Key4hep", "-u", "https://key4hep.github.io/key4hep-doc/", - "-x", ".github/*", ".pre-commit-config.yaml", "README.md", "doc/ReleaseNotes.md", "*.xml", "-f"] + "-x", ".github/*", ".pre-commit-config.yaml", "README.md", "doc/ReleaseNotes.md", "*.xml", "*.ipynb", "-f"] diff --git a/CMakeLists.txt b/CMakeLists.txt index 1949b6f..ee30546 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,6 +18,11 @@ limitations under the License. ]] CMAKE_MINIMUM_REQUIRED(VERSION 3.12) +# DD4hep 1.36 still uses CMake's FindBoost module. Keep that compatibility +# mode explicitly on CMake 3.30+ instead of emitting a CMP0167 developer warning. +if(POLICY CMP0167) + cmake_policy(SET CMP0167 OLD) +endif() project(SoftwareTutorials LANGUAGES CXX) diff --git a/DD4hepTutorials/CMakeLists.txt b/DD4hepTutorials/CMakeLists.txt index 16685f5..9e9cf29 100644 --- a/DD4hepTutorials/CMakeLists.txt +++ b/DD4hepTutorials/CMakeLists.txt @@ -17,6 +17,12 @@ See the License for the specific language governing permissions and limitations under the License. ]] CMAKE_MINIMUM_REQUIRED(VERSION 3.12) + +# DD4hep 1.36 still relies on CMake's legacy FindBoost compatibility mode. +if(POLICY CMP0167) + cmake_policy(SET CMP0167 OLD) +endif() + project(DD4hepTutorials) include(GNUInstallDirs) diff --git a/DD4hepTutorials/README.md b/DD4hepTutorials/README.md index 55696fe..ca443cd 100644 --- a/DD4hepTutorials/README.md +++ b/DD4hepTutorials/README.md @@ -18,17 +18,88 @@ limitations under the License. --> # DD4hepTutorials for DRDCalo -Welcome to the DD4hep Tutorials of DRDCalo repository! +Hands-on exercises for DD4hep, a detector description framework: -This repository contains hands-on exercises to help you learn DD4hep, a powerful detector description framework. The exercises are organized into separate folders: + - **simplecalo1**: the fundamentals, by building a simplified calorimeter. + - **simplecalo2**: a custom Geant4 sensitive detector on top of it. - - **simplecalo1**: Learn the fundamentals of DD4hep by building your first simplified calorimeter. - - **simplecalo2**: Build on simplecalo1 and explore how to set up a custom Geant4 Sensitive Detector. - -For instructions on **simplecalo1** and **simplecalo2** follow [this presentation](https://indico.cern.ch/event/1618975/sessions/635708/attachments/3252510/5805581/DRDCaloDD4hepTutorial_April2026.pdf). +Follow [this presentation](https://indico.cern.ch/event/1618975/sessions/635708/attachments/3252510/5805581/DRDCaloDD4hepTutorial_April2026.pdf) for both. ## Previous editions of this tutorial - April 2026, DRDCalo Collaboration Meeting [presentation](https://indico.cern.ch/event/1618975/sessions/635708/attachments/3252510/5805581/DRDCaloDD4hepTutorial_April2026.pdf) - April 2025, DRDCalo Collaboration Meeting [presentation](https://indico.ijclab.in2p3.fr/event/11400/sessions/5873/attachments/25413/37372/DRD6DD4hepTutorial_April2025.pdf) +## Analysing the output + +| What | Notebook | +|---|---| +| Section 1: cell energy sum, Gaussian fit, energy resolution | `notebooks/plot_cell_energy_sum.ipynb` | +| Hands-on 6: hits, layers, lateral shape, contributions | `notebooks/readEdm4hep.ipynb` | + +Hands-on 6 has six questions to complete; the finished version sits next to it as +`notebooks/readEdm4hepSolution.ipynb`. + +Each notebook prefers the full simulation — `simplecalo1.root` and `simplecalo2.root`, which the +steering files always write into this directory whatever directory `ddsim` was launched from — +and otherwise falls back to the bundled 10-event `data/simplecalo2_sample.root`, saying so when +it does. `SIMPLECALO1_FILE` and `SIMPLECALO2_FILE` override the choice; the rules live in +`python/drdcalo_tutorials/__init__.py`. + +Hands-on 6 reads the cell size and grid dimensions out of `simplecalo2/compact/simplecalo2.xml`, +the same constants `sc2_solution1.h` places the cells from, so changing `CellX` and re-running +`ddsim` needs no code edit. Its last cell checks a decoded cell index against the stored hit +position, which catches an input file produced before the geometry was changed. + +### Running the notebooks + +Open `SoftwareTutorials` itself as the VS Code workspace folder — not a parent directory — then +**Select Kernel** → **Python Environments** → `.venv/bin/python`. VS Code applies +`${workspaceFolder}/.env`, which `setup.sh` fills with the key4hep runtime. A `No module named +'awkward'` error means it did not, almost always because the workspace folder is a level too +high; see the main README. + +For JupyterLab instead, run `jupyter lab --no-browser --port 8888` from a shell that has sourced +`setup.sh`, and forward that port. + +## Viewing the geometry + +`geoWebDisplay simplecalo1/compact/simplecalo1.xml` starts a web server and asks the operating +system to open a browser. On a remote machine that launch usually fails silently — the server is +still running, on the port it printed: + +``` +Info in : Starting HTTP server on port 127.0.0.1:9427 +``` + +Leave the `root [0]` prompt open, since quitting it stops the server, and forward that port to +your own machine — either through the VS Code **PORTS** panel, or: + +``` bash +# on your laptop. Use the exact node name: lxplus.cern.ch is load balanced and would +# send you to a different machine, where nothing is listening. +ssh -N -L 9427:127.0.0.1:9427 @lxplus8sXX.cern.ch +``` + +Then browse to ****. The `/win1/` matters: ROOT serves each GUI panel +as a separate named window, so the bare URL returns *404*. `/win2/` is the volume hierarchy +browser, useful on its own during Hands-on 1 and 2. + +> ROOT also appends a single-use key to the URL and hands it only to the browser it launches +> itself, so a typed address returns *404* as well. To type one, put `WebGui.OnetimeKey: no` into +> a `.rootrc` in the directory you start `geoWebDisplay` from. + +`geoWebDisplay` also drops a `viewer.cxx` in the current directory; it is a ROOT by-product, +ignored by git, and can be deleted. + +### Without port forwarding + +Export the geometry and look at it locally: + +``` bash +geoConverter -compact2tgeo -input simplecalo1/compact/simplecalo1.xml -output simplecalo1_geo.root +# or -compact2gdml for a GDML file +``` + +Open the result in ROOT, or drag it onto . + Happy coding! :rocket: diff --git a/DD4hepTutorials/data/README.md b/DD4hepTutorials/data/README.md new file mode 100644 index 0000000..77b9821 --- /dev/null +++ b/DD4hepTutorials/data/README.md @@ -0,0 +1,59 @@ + +# Bundled EDM4hep sample + +`simplecalo2_sample.root` is a deterministic 10-event input shared by the DD4hep analysis tools +and the Gaudi exercises. It contains `EventHeader`, `MCParticles`, `simplecaloRO`, and +`simplecaloROContributions` and was generated with the completed simplecalo2 geometry using +key4hep release `2026-04-08` and random seed `12345`. + +All lengths are in **mm** and all times in **ns**, as EDM4hep specifies: `simplecaloRO.position` +is the centre of the cell, `simplecaloROContributions.stepPosition` the position of the +individual Geant4 step that contributed. Energies are in GeV for hits and contributions. + +The small sample lets the notebook run after a fresh clone. It is not intended to replace the +full 500-event simulation when producing statistically meaningful plots — the notebooks say so +out loud when they fall back to it. After completing the cell-placement exercise, regenerate a +full `DD4hepTutorials/simplecalo2.root` with: + +```bash +ddsim --steeringFile DD4hepTutorials/simplecalo2/sc2SteeringFile.py +``` + +The steering file names both its compact file and its output relative to itself, so the command +above works from any directory and always writes `DD4hepTutorials/simplecalo2.root`, which is +where the notebooks look. Pass `--outputFile` to put it somewhere else, and point the notebooks +at it with `SIMPLECALO2_FILE`. + +To reproduce the bundled fixture after enabling the completed cell-placement implementation: + +```bash +ddsim --steeringFile DD4hepTutorials/simplecalo2/sc2SteeringFile.py \ + --numberOfEvents 10 \ + --outputFile DD4hepTutorials/data/simplecalo2_sample.root \ + --random.seed 12345 \ + --random.enableEventSeed +``` + +The geometry it was produced with is the one in `simplecalo2/compact/simplecalo2.xml` at the time: +10 x 10 cells of 10 cm. Change `CellX` there and this file no longer matches — the last cell of +`readEdm4hepSolution.ipynb` compares a decoded cell index against the stored hit position and will +say so. + +SHA-256: `a4faa81b48a13e48f4363224e039c41f949d3452acdd223ebaba61b9f1b672b2` diff --git a/DD4hepTutorials/data/simplecalo2_sample.root b/DD4hepTutorials/data/simplecalo2_sample.root new file mode 100644 index 0000000..c833ab3 Binary files /dev/null and b/DD4hepTutorials/data/simplecalo2_sample.root differ diff --git a/DD4hepTutorials/notebooks/plot_cell_energy_sum.ipynb b/DD4hepTutorials/notebooks/plot_cell_energy_sum.ipynb new file mode 100644 index 0000000..85e845b --- /dev/null +++ b/DD4hepTutorials/notebooks/plot_cell_energy_sum.ipynb @@ -0,0 +1,304 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b4338f74", + "metadata": {}, + "source": [ + "\n", + "\n", + "# Section 1: cell energy sum and energy resolution\n", + "\n", + "Sums the energy of every calorimeter cell in each event, fits the peak with a Gaussian, and reports the sampling calorimeter's **energy resolution** σ/µ.\n", + "\n", + "Uproot reads the EDM4hep file, Awkward Array and NumPy do the arithmetic, SciPy fits, and Matplotlib draws.\n", + "\n", + "> **Hands-on 3** lives at the bottom of this notebook: change the sampling fraction and frequency in the compact file, re-run `ddsim`, and watch σ/µ move." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ac768d65", + "metadata": {}, + "outputs": [], + "source": [ + "import awkward as ak\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import uproot\n", + "from scipy.optimize import curve_fit\n", + "\n", + "from drdcalo_tutorials import simplecalo1_input\n", + "\n", + "# One blue for the data, one orange for the model, and three greys for everything\n", + "# that is not data. Keeping ink and marks apart is what stops a plot looking busy.\n", + "BLUE, ORANGE = \"#2a78d6\", \"#eb6834\"\n", + "INK, INK_SOFT, GRID, SURFACE = \"#0b0b0b\", \"#52514e\", \"#dcdcd8\", \"#fcfcfb\"\n", + "\n", + "plt.rcParams.update({\n", + " \"figure.facecolor\": SURFACE,\n", + " \"axes.facecolor\": SURFACE,\n", + " \"axes.edgecolor\": GRID,\n", + " \"axes.labelcolor\": INK_SOFT,\n", + " \"axes.titlecolor\": INK,\n", + " \"xtick.color\": INK_SOFT,\n", + " \"ytick.color\": INK_SOFT,\n", + " \"xtick.labelcolor\": INK_SOFT,\n", + " \"ytick.labelcolor\": INK_SOFT,\n", + " \"font.size\": 11,\n", + " \"figure.dpi\": 120,\n", + "})" + ] + }, + { + "cell_type": "markdown", + "id": "1d20aa4d", + "metadata": {}, + "source": [ + "## Open the event data\n", + "\n", + "`simplecalo1_input()` prefers the `simplecalo1.root` you produced with\n", + "\n", + "```bash\n", + "ddsim --steeringFile simplecalo1/sc1SteeringFile.py\n", + "```\n", + "\n", + "and falls back to the bundled 10-event sample so that this notebook also runs straight after a clone. `SIMPLECALO1_FILE` overrides both." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1d2a8061", + "metadata": {}, + "outputs": [], + "source": [ + "input_file = simplecalo1_input()\n", + "COLLECTION = \"simplecaloRO\"\n", + "\n", + "events = uproot.open(input_file, handler=uproot.MultithreadedFileSource)[\"events\"]\n", + "n_events = events.num_entries\n", + "print(f\"Reading {n_events} events from {input_file}\")\n", + "\n", + "# A Gaussian fitted to a handful of events converges and prints a number, and that\n", + "# number is noise. Say so rather than letting the 10-event fallback look like a\n", + "# measurement of the calorimeter.\n", + "MIN_EVENTS_FOR_A_RESOLUTION = 100\n", + "if n_events < MIN_EVENTS_FOR_A_RESOLUTION:\n", + " print(\n", + " f\"\\nWARNING: {n_events} events is far too few to quote an energy resolution.\\n\"\n", + " f\" Generate the full simplecalo1 simulation (see the cell above)\\n\"\n", + " f\" and re-run this notebook before believing the sigma/mu below.\"\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "6ffd9dce", + "metadata": {}, + "source": [ + "## Sum the cells, event by event\n", + "\n", + "Uproot hands back one nested list of cell energies per event. Summing on `axis=1` collapses each of those lists to a single number and keeps the event boundary, so no Python loop over events is needed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "49321c04", + "metadata": {}, + "outputs": [], + "source": [ + "cell_energies = events[f\"{COLLECTION}/{COLLECTION}.energy\"].array(library=\"ak\")\n", + "total_energy = ak.to_numpy(ak.sum(cell_energies, axis=1)).astype(float)\n", + "\n", + "if total_energy.size == 0:\n", + " raise SystemExit(f\"No events found in {input_file}\")\n", + "\n", + "print(f\"{ak.sum(ak.num(cell_energies))} cells in {n_events} events\")\n", + "print(f\"Mean energy sum: {total_energy.mean():.4f} GeV\")" + ] + }, + { + "cell_type": "markdown", + "id": "2151dd1f", + "metadata": {}, + "source": [ + "## Histogram and Gaussian fit\n", + "\n", + "The binning follows the data instead of a fixed axis: the energy scale depends on the sampling fraction, which Hands-on 3 invites you to change.\n", + "\n", + "Only the core of the peak is fitted — the highest bin ± 3 RMS — so that the tails, where a shower leaked out of the calorimeter, do not pull the width. Bin uncertainties are Poisson, which is what makes the reported σ uncertainty and the χ² meaningful." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dbed0b60", + "metadata": {}, + "outputs": [], + "source": [ + "BINS = 50\n", + "\n", + "\n", + "def gaussian(x, amplitude, mean, sigma):\n", + " return amplitude * np.exp(-0.5 * ((x - mean) / sigma) ** 2)\n", + "\n", + "\n", + "# The binning follows the data rather than a fixed 0-4 GeV axis: the energy scale\n", + "# depends on the sampling fraction, which Hands-on 3 invites you to change, and a\n", + "# range anchored at zero would waste most of the plot on an empty region.\n", + "margin = 0.05 * (total_energy.max() - total_energy.min()) or 1.0\n", + "counts, edges = np.histogram(\n", + " total_energy, bins=BINS, range=(total_energy.min() - margin, total_energy.max() + margin)\n", + ")\n", + "centres = 0.5 * (edges[:-1] + edges[1:])\n", + "\n", + "peak = centres[np.argmax(counts)]\n", + "half_window = 3 * total_energy.std()\n", + "\n", + "# Empty bins carry no Poisson weight, so they are left out of the fit rather than\n", + "# given an arbitrary uncertainty.\n", + "in_fit = (np.abs(centres - peak) <= half_window) & (counts > 0)\n", + "\n", + "popt, pcov = curve_fit(\n", + " gaussian,\n", + " centres[in_fit],\n", + " counts[in_fit],\n", + " p0=(counts.max(), peak, total_energy.std()),\n", + " sigma=np.sqrt(counts[in_fit]),\n", + " absolute_sigma=True,\n", + ")\n", + "\n", + "amplitude, mean, sigma = popt\n", + "sigma = abs(sigma) # the Gaussian is symmetric in sigma, so the fit may return either sign\n", + "mean_error, sigma_error = np.sqrt(np.diag(pcov))[1:]\n", + "\n", + "pulls = (counts[in_fit] - gaussian(centres[in_fit], *popt)) / np.sqrt(counts[in_fit])\n", + "chi2, ndf = float(np.sum(pulls**2)), int(in_fit.sum() - len(popt))\n", + "resolution = sigma / mean" + ] + }, + { + "cell_type": "markdown", + "id": "4e565495", + "metadata": {}, + "source": [ + "## Result" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "510d473c", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(6, 5))\n", + "\n", + "# The distribution: a soft fill for the shape, a crisp outline for the edges, and\n", + "# Poisson bars so the reader can see what the fit was actually weighted by\n", + "ax.stairs(counts, edges, fill=True, facecolor=BLUE, alpha=0.16, edgecolor=\"none\", zorder=1)\n", + "ax.stairs(counts, edges, color=BLUE, linewidth=1.6, label=\"Simulation\", zorder=2)\n", + "ax.errorbar(centres, counts, yerr=np.sqrt(counts), fmt=\"none\",\n", + " ecolor=BLUE, elinewidth=1.0, alpha=0.5, zorder=3)\n", + "\n", + "# Drawing the curve over the fit window only, so that where it stops is also\n", + "# where the fit stopped looking\n", + "fit_curve_x = np.linspace(peak - half_window, peak + half_window, 400)\n", + "ax.plot(fit_curve_x, gaussian(fit_curve_x, *popt), color=ORANGE, linewidth=2,\n", + " label=\"Gaussian fit\", zorder=4)\n", + "\n", + "ax.set_title(\"Total cell energy per event\", loc=\"left\", fontsize=13,\n", + " fontweight=\"semibold\", pad=26)\n", + "ax.text(0, 1.04, f\"simplecalo1 · {n_events} events\", transform=ax.transAxes,\n", + " fontsize=10, color=INK_SOFT)\n", + "\n", + "ax.text(0.98, 0.94,\n", + " \"\\n\".join([\n", + " rf\"$\\mu$ = {mean:.4f} $\\pm$ {mean_error:.4f} GeV\",\n", + " rf\"$\\sigma$ = {sigma:.4f} $\\pm$ {sigma_error:.4f} GeV\",\n", + " rf\"$\\sigma/\\mu$ = {100 * resolution:.2f} %\",\n", + " rf\"$\\chi^2$/ndf = {chi2:.1f} / {ndf}\",\n", + " ]),\n", + " transform=ax.transAxes, ha=\"right\", va=\"top\",\n", + " fontsize=10, color=INK_SOFT, linespacing=1.7)\n", + "\n", + "# Say it on the figure too, so that a screenshot cannot be mistaken for a result\n", + "if n_events < MIN_EVENTS_FOR_A_RESOLUTION:\n", + " ax.text(0.5, 0.55, f\"only {n_events} events\\nnot a measurement\",\n", + " transform=ax.transAxes, ha=\"center\", va=\"center\", rotation=18,\n", + " fontsize=20, fontweight=\"bold\", color=ORANGE, alpha=0.30, zorder=5)\n", + "\n", + "ax.set_xlabel(\"Cell energy sum [GeV]\")\n", + "ax.set_ylabel(\"Events\")\n", + "ax.set_xlim(edges[0], edges[-1])\n", + "# Headroom, so that the legend and the fit results never sit on top of the peak\n", + "ax.set_ylim(0, 1.4 * counts.max())\n", + "ax.set_axisbelow(True)\n", + "ax.spines[\"top\"].set_visible(False)\n", + "ax.spines[\"right\"].set_visible(False)\n", + "ax.legend(frameon=False, loc=\"upper left\", fontsize=10, labelcolor=INK_SOFT)\n", + "\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a778df97", + "metadata": {}, + "outputs": [], + "source": [ + "print(f\"Fitted mean {mean:.4f} +/- {mean_error:.4f} GeV\")\n", + "print(f\"Fitted sigma {sigma:.4f} +/- {sigma_error:.4f} GeV\")\n", + "print(f\"chi2 / ndf {chi2:.1f} / {ndf}\")\n", + "print()\n", + "print(f\"Energy resolution sigma/mean = {100 * resolution:.2f} %\")\n", + "\n", + "if n_events < MIN_EVENTS_FOR_A_RESOLUTION:\n", + " print(f\" ... from {n_events} events, so this is noise, not a resolution\")" + ] + }, + { + "cell_type": "markdown", + "id": "84087c67", + "metadata": {}, + "source": [ + "## Hands-on 3 — sampling fraction and sampling frequency\n", + "\n", + "The number above is the whole point of the exercise. Change the calorimeter in `simplecalo1/compact/simplecalo1.xml` and see how it responds:\n", + "\n", + "- **Sampling fraction** — the share of a layer that is active. Move `AbsLayerZ` and `SensLayerZ` while keeping `CaloLayerZ = AbsLayerZ + SensLayerZ`. A thicker active layer collects more of the shower, so the visible energy grows.\n", + "- **Sampling frequency** — how finely the shower is sampled. Shrink `CaloLayerZ` and raise `LayersNumber` so the calorimeter keeps its total depth. More, thinner layers sample the shower more often.\n", + "\n", + "Then re-run the simulation and this notebook:\n", + "\n", + "```bash\n", + "ddsim --steeringFile simplecalo1/sc1SteeringFile.py\n", + "```\n", + "\n", + "Sampling calorimeters are dominated by sampling fluctuations, so σ/µ should improve roughly as the square root of the sampling frequency. Watch the fitted µ as well: it tells you how much of the 50 GeV actually ends up in the active material.\n", + "\n", + "> Nothing is hard-coded to the default geometry — the histogram range follows the data — so the fit keeps working as the energy scale moves." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/DD4hepTutorials/notebooks/readEdm4hep.ipynb b/DD4hepTutorials/notebooks/readEdm4hep.ipynb new file mode 100644 index 0000000..14d4498 --- /dev/null +++ b/DD4hepTutorials/notebooks/readEdm4hep.ipynb @@ -0,0 +1,329 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "ba1a753f", + "metadata": {}, + "source": [ + "\n", + "\n", + "# Hands-on 6: reading EDM4hep hits and contributions\n", + "\n", + "This exercise reads the **simplecalo2** EDM4hep file in Python and builds the total energy spectrum, longitudinal shower profile, lateral shower shape, and contribution timing distribution.\n", + "\n", + "There are six questions, marked `Q1` to `Q6`. The completed notebook is `readEdm4hepSolution.ipynb`. Uproot reads the events, DD4hep decodes cell IDs, and Awkward Array, NumPy, and Matplotlib handle the analysis and plots." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "263ba7d5", + "metadata": {}, + "outputs": [], + "source": [ + "import awkward as ak\n", + "import dd4hep\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.colors import LogNorm\n", + "import numpy as np\n", + "import uproot\n", + "\n", + "from drdcalo_tutorials import SIMPLECALO2_COMPACT, compact_constants, simplecalo2_input" + ] + }, + { + "cell_type": "markdown", + "id": "7a18acd1", + "metadata": {}, + "source": [ + "## Open the event data\n", + "\n", + "The generated 500-event file is used when available; otherwise the notebook uses the bundled 10-event sample. Set `SIMPLECALO2_FILE` to choose another compatible file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f9ad9ae3", + "metadata": {}, + "outputs": [], + "source": [ + "input_file = simplecalo2_input()\n", + "\n", + "events = uproot.open(input_file, handler=uproot.MultithreadedFileSource)[\"events\"]\n", + "print(f\"Reading {events.num_entries} events from {input_file}\")" + ] + }, + { + "cell_type": "markdown", + "id": "cc61ae39", + "metadata": {}, + "source": [ + "## Q1 — decode the cell ID\n", + "\n", + "Each hit carries a 64-bit `cellID`. Which encoding string belongs to `simplecaloRO`? Find the `` element in the `` block of `simplecalo2/compact/simplecalo2.xml`, then pass it to DD4hep's native bit-field decoder." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "931839c9", + "metadata": {}, + "outputs": [], + "source": [ + "# Q1\n", + "ENCODING = \"FILL ME\"\n", + "decoder = dd4hep.core.DDSegmentation.BitFieldCoder(ENCODING)" + ] + }, + { + "cell_type": "markdown", + "id": "f0ceffab", + "metadata": {}, + "source": [ + "## Q2 and Q3 — read and decode the hits\n", + "\n", + "Uproot exposes each EDM4hep member as an Awkward array with one nested list per event. For Q2, replace `FILL ME` with the readout collection name. For Q3, use the `decode` helper below to extract `calolayer`, `abslayer`, and `cellid` from the flattened IDs." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "243e6524", + "metadata": {}, + "outputs": [], + "source": [ + "# Q2\n", + "COLLECTION = \"FILL ME\"\n", + "\n", + "def read(field, collection=COLLECTION):\n", + " return events[f\"{collection}/{collection}.{field}\"].array(library=\"ak\")\n", + "\n", + "cell_ids = read(\"cellID\")\n", + "hit_energies = read(\"energy\")\n", + "hit_x = read(\"position.x\")\n", + "hit_y = read(\"position.y\")\n", + "\n", + "# Q3\n", + "# decoder.get(cellID, \"field\") is DD4hep's own decoder, the same call the\n", + "# cross-check at the bottom of this notebook makes. int() because uproot hands\n", + "# back NumPy integers and DD4hep wants a plain Python one.\n", + "flat_cell_ids = ak.to_numpy(ak.flatten(cell_ids))\n", + "\n", + "def decode(field):\n", + " return np.array([decoder.get(int(cell_id), field) for cell_id in flat_cell_ids])\n", + "\n", + "calo_layer = ... # FILL ME: decode(\"...\")\n", + "abs_layer = ... # FILL ME\n", + "sub_cell_id = ... # FILL ME\n", + "\n", + "# The geometry is read back from the same compact file the simulation used, so\n", + "# that changing CellX or SensLayerX in simplecalo2.xml changes this analysis with\n", + "# it, exactly as it changes the cell placement in sc2_solution1.h.\n", + "geometry = compact_constants(SIMPLECALO2_COMPACT)\n", + "N_LAYERS = int(geometry[\"LayersNumber\"])\n", + "CELL_X, CELL_Y = geometry[\"CellX\"], geometry[\"CellY\"] # mm\n", + "N_CELLS_X = round(geometry[\"SensLayerX\"] / CELL_X)\n", + "N_CELLS_Y = round(geometry[\"SensLayerY\"] / CELL_Y)\n", + "print(f\"{N_CELLS_X} x {N_CELLS_Y} cells of {CELL_X:.0f} mm x {CELL_Y:.0f} mm in {N_LAYERS} layers\")\n", + "\n", + "# sc2_solution1.h runs the outer loop over y and the inner one over x, so\n", + "# cellid = N_CELLS_X * iY + iX. That is the convention undone here.\n", + "cell_x = sub_cell_id % N_CELLS_X\n", + "cell_y = sub_cell_id // N_CELLS_X" + ] + }, + { + "cell_type": "markdown", + "id": "277dadb5", + "metadata": {}, + "source": [ + "## Q4 — total energy and longitudinal profile\n", + "\n", + "Use `ak.sum(..., axis=1)` to keep one total per event. Flatten the hit energies with `ak.flatten`, then use the decoded layer array and `np.histogram(..., weights=...)` for the layer profile. Divide the profile by the number of events." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "91049fe8", + "metadata": {}, + "outputs": [], + "source": [ + "# Q4\n", + "n_events = events.num_entries\n", + "total_energies = ... # FILL ME\n", + "energy = ... # FILL ME\n", + "layer = calo_layer\n", + "\n", + "energy_max = 1.2 * total_energies.max() if total_energies.max() > 0 else 1.0\n", + "energy_counts, energy_edges = np.histogram(total_energies, bins=100, range=(0, energy_max))\n", + "layer_edges = np.arange(0.5, N_LAYERS + 1.5)\n", + "average_layer_energy, _ = ... # FILL ME\n", + "average_layer_energy /= n_events\n", + "\n", + "print(f\"Mean total energy: {total_energies.mean():.3f} GeV over {n_events} events\")" + ] + }, + { + "cell_type": "markdown", + "id": "2cda676f", + "metadata": {}, + "source": [ + "## Q5 — lateral shower shape\n", + "\n", + "The cell placement uses the outer loop for *y* and the inner loop for *x*, so `cellid = N_CELLS_X × iY + iX`. Build one weighted `np.histogram2d` per layer and divide it by the number of events." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dd9eaa0b", + "metadata": {}, + "outputs": [], + "source": [ + "# Q5\n", + "ix = cell_x\n", + "iy = cell_y\n", + "lateral_energy = []\n", + "\n", + "for layer_number in range(1, N_LAYERS + 1):\n", + " in_layer = layer == layer_number\n", + " histogram, _, _ = ... # FILL ME\n", + " lateral_energy.append(histogram / n_events)\n", + "\n", + "lateral_energy = np.stack(lateral_energy)" + ] + }, + { + "cell_type": "markdown", + "id": "fd3fd6ed", + "metadata": {}, + "source": [ + "## Q6 — contribution timing\n", + "\n", + "EDM4hep stores the hit contributions in the companion collection `Contributions`. Read its `time` and `energy` fields, flatten them, and make a 50-bin weighted histogram from 0 to 10 ns. Each contribution also carries `stepPosition.x/y/z` and `stepLength`, the position and length of the Geant4 step that produced it, in mm." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ed8f19f", + "metadata": {}, + "outputs": [], + "source": [ + "# Q6\n", + "CONTRIBUTIONS = f\"{COLLECTION}Contributions\"\n", + "contribution_times = ... # FILL ME\n", + "contribution_energies = ... # FILL ME\n", + "time_energy, time_edges = ... # FILL ME\n", + "print(f\"Read {len(contribution_times)} contributions\")" + ] + }, + { + "cell_type": "markdown", + "id": "e5062c97", + "metadata": {}, + "source": [ + "## Results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "87950c5d", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", + "axes[0].stairs(energy_counts, energy_edges, fill=True)\n", + "axes[0].set(xlabel=\"Energy [GeV]\", ylabel=\"Events\", title=\"Total energy deposit\")\n", + "axes[1].bar(range(1, N_LAYERS + 1), average_layer_energy)\n", + "axes[1].set(xlabel=\"Layer\", ylabel=\"Energy [GeV]\", title=\"Average energy deposit per layer\")\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cb704722", + "metadata": {}, + "outputs": [], + "source": [ + "busiest = int(np.argmax(lateral_energy.sum(axis=(1, 2))))\n", + "lateral = lateral_energy[busiest]\n", + "print(f\"Most energetic layer: {busiest + 1}\")\n", + "print(f\"Hottest cell: {100 * lateral.max() / lateral.sum():.1f}% of that layer's energy\")\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", + "# Brass has a Moliere radius of about 2 cm, small against cells this size, so the\n", + "# impact cell takes almost everything. A log scale is what makes the tails visible.\n", + "image = axes[0].imshow(\n", + " np.ma.masked_less_equal(lateral.T, 0), origin=\"lower\",\n", + " extent=(0, N_CELLS_X, 0, N_CELLS_Y), aspect=\"equal\",\n", + " norm=LogNorm(vmin=lateral.max() * 1e-5, vmax=lateral.max()),\n", + ")\n", + "axes[0].set(xlabel=r\"$i_X$\", ylabel=r\"$i_Y$\", title=f\"Lateral shape, layer {busiest + 1}\")\n", + "fig.colorbar(image, ax=axes[0], label=\"Average energy [GeV]\")\n", + "axes[1].stairs(time_energy, time_edges, fill=True)\n", + "axes[1].set(xlabel=\"Time [ns]\", ylabel=\"Energy [GeV]\", title=\"Hit contribution timing\")\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "10fd19ea", + "metadata": {}, + "source": [ + "## Cross-check the cell indices\n", + "\n", + "The hottest hit in the first event should have the same stored position as the cell centre implied by its decoded indices." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "67a118e2", + "metadata": {}, + "outputs": [], + "source": [ + "hottest = int(np.argmax(ak.to_numpy(hit_energies[0])))\n", + "sub = decoder.get(int(cell_ids[0][hottest]), \"cellid\")\n", + "i_x, i_y = sub % N_CELLS_X, sub // N_CELLS_X\n", + "\n", + "# The same arithmetic as sc2_solution1.h. EDM4hep stores positions in mm and\n", + "# compact_constants() returns the XML dimensions in mm, so the two compare directly.\n", + "implied_x = -geometry[\"SensLayerX\"] / 2 + CELL_X / 2 + i_x * CELL_X\n", + "implied_y = geometry[\"SensLayerY\"] / 2 - CELL_Y / 2 - i_y * CELL_Y\n", + "stored_x, stored_y = float(hit_x[0][hottest]), float(hit_y[0][hottest])\n", + "\n", + "print(f\"hottest cell: cellid={sub} -> iX={i_x}, iY={i_y}\")\n", + "print(f\" implied centre : x={implied_x:+7.1f} mm, y={implied_y:+7.1f} mm\")\n", + "print(f\" stored position: x={stored_x:+7.1f} mm, y={stored_y:+7.1f} mm\")\n", + "\n", + "if np.isclose(implied_x, stored_x) and np.isclose(implied_y, stored_y):\n", + " print(\" -> match: the cell indices are decoded correctly\")\n", + "else:\n", + " print(\" -> MISMATCH: the file being read was not produced by the geometry now\")\n", + " print(\" in simplecalo2.xml. Either it predates Hands-on 4, or the compact\")\n", + " print(\" file was edited afterwards and ddsim has not been re-run since.\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/DD4hepTutorials/notebooks/readEdm4hepSolution.ipynb b/DD4hepTutorials/notebooks/readEdm4hepSolution.ipynb new file mode 100644 index 0000000..ccfc798 --- /dev/null +++ b/DD4hepTutorials/notebooks/readEdm4hepSolution.ipynb @@ -0,0 +1,338 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b5b02277", + "metadata": {}, + "source": [ + "\n", + "\n", + "# Hands-on 6: reading EDM4hep hits and contributions\n", + "\n", + "This is the completed version of `readEdm4hep.ipynb`. It reads the **simplecalo2** EDM4hep file in Python and builds the total energy spectrum, longitudinal shower profile, lateral shower shape, and contribution timing distribution.\n", + "\n", + "Uproot reads the events, DD4hep decodes cell IDs, and Awkward Array, NumPy, and Matplotlib handle the analysis and plots. No ROOT histograms, canvases, or output files are used." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2fa23983", + "metadata": {}, + "outputs": [], + "source": [ + "import awkward as ak\n", + "import dd4hep\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.colors import LogNorm\n", + "import numpy as np\n", + "import uproot\n", + "\n", + "from drdcalo_tutorials import SIMPLECALO2_COMPACT, compact_constants, simplecalo2_input" + ] + }, + { + "cell_type": "markdown", + "id": "47f07c0e", + "metadata": {}, + "source": [ + "## Open the event data\n", + "\n", + "The generated 500-event file is used when available; otherwise the notebook uses the bundled 10-event sample. Set `SIMPLECALO2_FILE` to choose another compatible file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d4dea330", + "metadata": {}, + "outputs": [], + "source": [ + "input_file = simplecalo2_input()\n", + "\n", + "events = uproot.open(input_file, handler=uproot.MultithreadedFileSource)[\"events\"]\n", + "print(f\"Reading {events.num_entries} events from {input_file}\")" + ] + }, + { + "cell_type": "markdown", + "id": "51f16f26", + "metadata": {}, + "source": [ + "## Q1 — decode the cell ID\n", + "\n", + "Each hit carries a 64-bit `cellID`. The encoding string in `simplecalo2/compact/simplecalo2.xml` gives each field's width, from least to most significant bits. DD4hep's native bit-field decoder interprets it for us." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c82b1574", + "metadata": {}, + "outputs": [], + "source": [ + "# Q1: the element of the block in simplecalo2.xml\n", + "ENCODING = \"calolayer:5,abslayer:1,cellid:10\"\n", + "decoder = dd4hep.core.DDSegmentation.BitFieldCoder(ENCODING)\n", + "print(\"fields:\", [field.name() for field in decoder.fields()])" + ] + }, + { + "cell_type": "markdown", + "id": "e8e4cd52", + "metadata": {}, + "source": [ + "## Q2 and Q3 — read and decode the hits\n", + "\n", + "Uproot exposes each EDM4hep member as an Awkward array with one nested list per event. The hit collection is named after the detector readout. The `decode` helper below asks DD4hep's decoder for one field of every hit." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c98f8199", + "metadata": {}, + "outputs": [], + "source": [ + "# Q2: the collection is named after the readout in the compact file\n", + "COLLECTION = \"simplecaloRO\"\n", + "\n", + "def read(field, collection=COLLECTION):\n", + " return events[f\"{collection}/{collection}.{field}\"].array(library=\"ak\")\n", + "\n", + "cell_ids = read(\"cellID\")\n", + "hit_energies = read(\"energy\")\n", + "hit_x = read(\"position.x\")\n", + "hit_y = read(\"position.y\")\n", + "\n", + "# Q3: decode the flattened cell IDs\n", + "# decoder.get(cellID, \"field\") is DD4hep's own decoder, the same call the\n", + "# cross-check at the bottom of this notebook makes. int() because uproot hands\n", + "# back NumPy integers and DD4hep wants a plain Python one.\n", + "flat_cell_ids = ak.to_numpy(ak.flatten(cell_ids))\n", + "\n", + "def decode(field):\n", + " return np.array([decoder.get(int(cell_id), field) for cell_id in flat_cell_ids])\n", + "\n", + "calo_layer = decode(\"calolayer\")\n", + "abs_layer = decode(\"abslayer\")\n", + "sub_cell_id = decode(\"cellid\")\n", + "\n", + "# The geometry is read back from the same compact file the simulation used, so\n", + "# that changing CellX or SensLayerX in simplecalo2.xml changes this analysis with\n", + "# it, exactly as it changes the cell placement in sc2_solution1.h.\n", + "geometry = compact_constants(SIMPLECALO2_COMPACT)\n", + "N_LAYERS = int(geometry[\"LayersNumber\"])\n", + "CELL_X, CELL_Y = geometry[\"CellX\"], geometry[\"CellY\"] # mm\n", + "N_CELLS_X = round(geometry[\"SensLayerX\"] / CELL_X)\n", + "N_CELLS_Y = round(geometry[\"SensLayerY\"] / CELL_Y)\n", + "print(f\"{N_CELLS_X} x {N_CELLS_Y} cells of {CELL_X:.0f} mm x {CELL_Y:.0f} mm in {N_LAYERS} layers\")\n", + "\n", + "# sc2_solution1.h runs the outer loop over y and the inner one over x, so\n", + "# cellid = N_CELLS_X * iY + iX. That is the convention undone here.\n", + "cell_x = sub_cell_id % N_CELLS_X\n", + "cell_y = sub_cell_id // N_CELLS_X\n", + "print(f\"Read {ak.sum(ak.num(hit_energies))} hits\")" + ] + }, + { + "cell_type": "markdown", + "id": "8359d3aa", + "metadata": {}, + "source": [ + "## Q4 — total energy and longitudinal profile\n", + "\n", + "Summing on `axis=1` keeps the event boundary and gives one total per event. Flattening removes that boundary so the hits can be histogrammed by layer. As in the original analysis, the layer histogram is divided by the number of events." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6771c86c", + "metadata": {}, + "outputs": [], + "source": [ + "# Q4\n", + "n_events = events.num_entries\n", + "total_energies = ak.to_numpy(ak.sum(hit_energies, axis=1))\n", + "energy = ak.to_numpy(ak.flatten(hit_energies)).astype(float)\n", + "layer = calo_layer\n", + "\n", + "energy_max = 1.2 * total_energies.max() if total_energies.max() > 0 else 1.0\n", + "energy_counts, energy_edges = np.histogram(total_energies, bins=100, range=(0, energy_max))\n", + "layer_edges = np.arange(0.5, N_LAYERS + 1.5)\n", + "average_layer_energy, _ = np.histogram(layer, bins=layer_edges, weights=energy)\n", + "average_layer_energy /= n_events\n", + "\n", + "print(f\"Mean total energy: {total_energies.mean():.3f} GeV over {n_events} events\")" + ] + }, + { + "cell_type": "markdown", + "id": "13080695", + "metadata": {}, + "source": [ + "## Q5 — lateral shower shape\n", + "\n", + "The cell placement uses the outer loop for *y* and the inner loop for *x*, so `cellid = N_CELLS_X × iY + iX`. Each slice below is the average energy map for one layer." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7d3e2054", + "metadata": {}, + "outputs": [], + "source": [ + "# Q5\n", + "ix = cell_x\n", + "iy = cell_y\n", + "lateral_energy = []\n", + "\n", + "for layer_number in range(1, N_LAYERS + 1):\n", + " in_layer = layer == layer_number\n", + " histogram, _, _ = np.histogram2d(\n", + " ix[in_layer], iy[in_layer],\n", + " bins=(N_CELLS_X, N_CELLS_Y), range=((0, N_CELLS_X), (0, N_CELLS_Y)),\n", + " weights=energy[in_layer],\n", + " )\n", + " lateral_energy.append(histogram / n_events)\n", + "\n", + "lateral_energy = np.stack(lateral_energy)" + ] + }, + { + "cell_type": "markdown", + "id": "d14545d2", + "metadata": {}, + "source": [ + "## Q6 — contribution timing\n", + "\n", + "EDM4hep stores the hit contributions in a companion collection. Reading that collection directly visits the same contribution objects as looping over `hit.getContributions()`, while keeping the analysis vectorized. Each contribution also carries `stepPosition.x/y/z` and `stepLength`, the position and length of the Geant4 step that produced it, in mm." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "77df6b8c", + "metadata": {}, + "outputs": [], + "source": [ + "# Q6\n", + "CONTRIBUTIONS = f\"{COLLECTION}Contributions\"\n", + "contribution_times = ak.to_numpy(ak.flatten(read(\"time\", CONTRIBUTIONS)))\n", + "contribution_energies = ak.to_numpy(ak.flatten(read(\"energy\", CONTRIBUTIONS))).astype(float)\n", + "time_energy, time_edges = np.histogram(\n", + " contribution_times, bins=50, range=(0, 10), weights=contribution_energies\n", + ")\n", + "print(f\"Read {len(contribution_times)} contributions\")" + ] + }, + { + "cell_type": "markdown", + "id": "f10b7aeb", + "metadata": {}, + "source": [ + "## Results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "02e57a5e", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", + "axes[0].stairs(energy_counts, energy_edges, fill=True)\n", + "axes[0].set(xlabel=\"Energy [GeV]\", ylabel=\"Events\", title=\"Total energy deposit\")\n", + "axes[1].bar(range(1, N_LAYERS + 1), average_layer_energy)\n", + "axes[1].set(xlabel=\"Layer\", ylabel=\"Energy [GeV]\", title=\"Average energy deposit per layer\")\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3825c360", + "metadata": {}, + "outputs": [], + "source": [ + "busiest = int(np.argmax(lateral_energy.sum(axis=(1, 2))))\n", + "lateral = lateral_energy[busiest]\n", + "print(f\"Most energetic layer: {busiest + 1}\")\n", + "print(f\"Hottest cell: {100 * lateral.max() / lateral.sum():.1f}% of that layer's energy\")\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n", + "# Brass has a Moliere radius of about 2 cm, small against cells this size, so the\n", + "# impact cell takes almost everything. A log scale is what makes the tails visible.\n", + "image = axes[0].imshow(\n", + " np.ma.masked_less_equal(lateral.T, 0), origin=\"lower\",\n", + " extent=(0, N_CELLS_X, 0, N_CELLS_Y), aspect=\"equal\",\n", + " norm=LogNorm(vmin=lateral.max() * 1e-5, vmax=lateral.max()),\n", + ")\n", + "axes[0].set(xlabel=r\"$i_X$\", ylabel=r\"$i_Y$\", title=f\"Lateral shape, layer {busiest + 1}\")\n", + "fig.colorbar(image, ax=axes[0], label=\"Average energy [GeV]\")\n", + "axes[1].stairs(time_energy, time_edges, fill=True)\n", + "axes[1].set(xlabel=\"Time [ns]\", ylabel=\"Energy [GeV]\", title=\"Hit contribution timing\")\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "f652b68b", + "metadata": {}, + "source": [ + "## Cross-check the cell indices\n", + "\n", + "The hottest hit in the first event should have the same stored position as the cell centre implied by its decoded indices." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7089eb95", + "metadata": {}, + "outputs": [], + "source": [ + "hottest = int(np.argmax(ak.to_numpy(hit_energies[0])))\n", + "sub = decoder.get(int(cell_ids[0][hottest]), \"cellid\")\n", + "i_x, i_y = sub % N_CELLS_X, sub // N_CELLS_X\n", + "\n", + "# The same arithmetic as sc2_solution1.h. EDM4hep stores positions in mm and\n", + "# compact_constants() returns the XML dimensions in mm, so the two compare directly.\n", + "implied_x = -geometry[\"SensLayerX\"] / 2 + CELL_X / 2 + i_x * CELL_X\n", + "implied_y = geometry[\"SensLayerY\"] / 2 - CELL_Y / 2 - i_y * CELL_Y\n", + "stored_x, stored_y = float(hit_x[0][hottest]), float(hit_y[0][hottest])\n", + "\n", + "print(f\"hottest cell: cellid={sub} -> iX={i_x}, iY={i_y}\")\n", + "print(f\" implied centre : x={implied_x:+7.1f} mm, y={implied_y:+7.1f} mm\")\n", + "print(f\" stored position: x={stored_x:+7.1f} mm, y={stored_y:+7.1f} mm\")\n", + "\n", + "if np.isclose(implied_x, stored_x) and np.isclose(implied_y, stored_y):\n", + " print(\" -> match: the cell indices are decoded correctly\")\n", + "else:\n", + " print(\" -> MISMATCH: the file being read was not produced by the geometry now\")\n", + " print(\" in simplecalo2.xml. Either it predates Hands-on 4, or the compact\")\n", + " print(\" file was edited afterwards and ddsim has not been re-run since.\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "DRDCalo Tutorial (3.13.8)", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.13.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/DD4hepTutorials/scripts/plot_cell_energy_sum.C b/DD4hepTutorials/scripts/plot_cell_energy_sum.C deleted file mode 100644 index 5d30964..0000000 --- a/DD4hepTutorials/scripts/plot_cell_energy_sum.C +++ /dev/null @@ -1,47 +0,0 @@ -// usage -// root -b -// .L plot_cell_energy_sum.C -// plot_cell_energy_sum(path_to_file) -// NB: in this case, a simple ROOT command in the interpreter (events->Draw("Sum$(simplecaloRO.energy)")) would give you the same information but the goal here is to teach how to properly analyze podio files (needed for more advanced analyses) - -#include - -#include "podio/Reader.h" - -#include "edm4hep/SimCalorimeterHitCollection.h" - -int plot_cell_energy_sum(std::string input_file) { - - // Use the podio reader to open the file (it could also be opened with plain ROOT but that would not preserve the podio functionalities) - auto reader = podio::makeReader(input_file); - - // Prepare a ROOT TH1 to store the cell energy sum - TH1F* cell_energy_sum_th1 = new TH1F("cell_energy_sum_th1", ";Cell energy sum [GeV]; Number of events", 100, 0, 4.0); - - // Fill the TH1 with the cell energy sum - for (size_t i = 0; i < reader.getEvents(); ++i) { - auto event = reader.readNextEvent(); - auto& calo_cells = event.get("simplecaloRO"); - float total_energy = 0.; - for (const auto& calo_cell : calo_cells) { - total_energy += calo_cell.getEnergy(); - } - cell_energy_sum_th1->Fill(total_energy); - } - - // Fit and draw the histogram - TCanvas* cell_energy_sum_canvas = new TCanvas("cell_energy_sum_canvas"); - float fit_range_min = cell_energy_sum_th1->GetXaxis()->GetBinCenter(cell_energy_sum_th1->GetMaximumBin()) - 3 * cell_energy_sum_th1->GetRMS(); - float fit_range_max = cell_energy_sum_th1->GetXaxis()->GetBinCenter(cell_energy_sum_th1->GetMaximumBin()) + 3 * cell_energy_sum_th1->GetRMS(); - cell_energy_sum_th1->Fit("gaus", "SQ", "", fit_range_min, fit_range_max); - gStyle->SetOptFit(1111); // to get the fit result displayed - cell_energy_sum_th1->Draw(); - // Print the canvas in png with a name based on the input file name - size_t pos = input_file.find(".root"); - if (pos != std::string::npos) { - input_file.replace(pos, 5, "_cell_energy_sum.png"); - } - cell_energy_sum_canvas->Print(input_file.c_str()); - - return 0; -} diff --git a/DD4hepTutorials/scripts/plot_cell_energy_sum.py b/DD4hepTutorials/scripts/plot_cell_energy_sum.py deleted file mode 100644 index cbcf2eb..0000000 --- a/DD4hepTutorials/scripts/plot_cell_energy_sum.py +++ /dev/null @@ -1,53 +0,0 @@ -# -# Copyright (c) 2020-2024 Key4hep-Project. -# -# This file is part of Key4hep. -# See https://key4hep.github.io/key4hep-doc/ for further info. -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -import argparse -import ROOT -# prevent ROOT to display things on the flight -ROOT.gROOT.SetBatch(ROOT.kTRUE) -from podio import root_io - -parser = argparse.ArgumentParser() -parser.add_argument("input_file", type=str, help="Path to the input file") -args = parser.parse_args() - -# Use the podio reader to open the file (it could also be opened with plain ROOT but that would not preserve the podio functionalities) -reader = root_io.Reader(args.input_file) - -# Prepare a ROOT TH1 to store the cell energy sum -cell_energy_sum_th1 = ROOT.TH1F("cell_energy_sum_th1", ";Cell energy sum [GeV]; Number of events", 100, 0, 4.0) - -# Fill the TH1 with the cell energy sum -for event in reader.get("events"): - calo_cells = event.get("simplecaloRO") - total_energy = sum([calo_cell.getEnergy() for calo_cell in calo_cells]) - cell_energy_sum_th1.Fill(total_energy) - -# Fit and draw the histogram -cell_energy_sum_canvas = ROOT.TCanvas("cell_energy_sum_canvas") -fit_range_min = cell_energy_sum_th1.GetXaxis().GetBinCenter(cell_energy_sum_th1.GetMaximumBin()) - 3 * cell_energy_sum_th1.GetRMS() -fit_range_max = cell_energy_sum_th1.GetXaxis().GetBinCenter(cell_energy_sum_th1.GetMaximumBin()) + 3 * cell_energy_sum_th1.GetRMS() -cell_energy_sum_th1.Fit("gaus", "SQ", "", fit_range_min, fit_range_max) -ROOT.gStyle.SetOptFit(1111) # to get the fit result displayed -cell_energy_sum_th1.Draw() -cell_energy_sum_canvas.Print(args.input_file.replace(".root", "_cell_energy_sum.png")) - - - - - diff --git a/DD4hepTutorials/scripts/readEdm4hep.cc b/DD4hepTutorials/scripts/readEdm4hep.cc deleted file mode 100644 index a6b8ed1..0000000 --- a/DD4hepTutorials/scripts/readEdm4hep.cc +++ /dev/null @@ -1,166 +0,0 @@ -/* - * Copyright (c) 2020-2024 Key4hep-Project. - * - * This file is part of Key4hep. - * See https://key4hep.github.io/key4hep-doc/ for further info. - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -// This is Hands-on 6, a more complex analysis for simplecalo2 -// Follow the code and complete the six questions along the code. -// Solutions and at the end of this file. - -// load required libraries for ROOT interpreter -R__LOAD_LIBRARY(libDDCore.so) -R__LOAD_LIBRARY(libpodio.so) -R__LOAD_LIBRARY(libpodioRootIO.so) -R__LOAD_LIBRARY(libedm4hep.so) - -#include "DDSegmentation/BitFieldCoder.h" - -// Declare the namespace for ROOT interpreter -namespace dd4hep { -namespace DDSegmentation { - class BitFieldCoder; -} -} // namespace dd4hep - -#include "podio/Frame.h" -#include "podio/ROOTReader.h" - -#include "edm4hep/SimCalorimeterHitCollection.h" - -#include "TFile.h" -#include "TH1F.h" -#include "TH2F.h" -#include "TROOT.h" - -// hands-on example to read edm4hep SimCalorimeterHitCollection and perform some basic analysis -void readEdm4hep(const std::string& filename) { - // Open the ROOT file using podio - std::vector files = {filename}; - auto reader = podio::ROOTReader(); - reader.openFiles(files); - - // event tree is always named "events" - unsigned int entries = reader.getEntries("events"); - - // bit field decoder - // Q1: which string should we use to decode the cellID for SimCalorimeterHitCollection in this file? - // (hint: check the cellID encoding in the DD4hep XML file used for simulation) - auto decoder = dd4hep::DDSegmentation::BitFieldCoder("/*fill me*/"); - - // define histograms - // add more histograms as needed for your analysis - TH1::SetDefaultSumw2(); // enable Sumw2 for all histograms - auto* hEnergy = new TH1F("hEnergy", "Total energy deposit;Energy (GeV);Events", 100, 0, 5.); // adjust range as needed - auto* hLayerEnergySum = - new TH1F("hLayerEnergySum", "Avg energy deposit for each active Layer;i-th layer;Energy (GeV)", 20, 0, 20); - auto* hContributionTime = new TH1F("hContributionTime", "Timing distribution of hit contributions;Time (ns);Entries", - 50, 0, 10); // adjust range as needed - - // lateral shower shape per layer - std::vector hLateralShapePerLayer; - for (int i = 1; i <= 20; i++) { // layer index starts from 1 - hLateralShapePerLayer.push_back(new TH2F(Form("hLateralShape_Layer%d", i), - Form("Lateral shower shape for Layer %d;iX;iY", i), 10, 0, 10, 10, 0, 10)); - } - - // loop over events - for (unsigned int iEvt = 0; iEvt < entries; iEvt++) { - if (iEvt % 100 == 0) - printf("Analyzing %dth event ...\n", iEvt); - - // read podio::Frame for the current event - auto aframe = podio::Frame(reader.readEntry("events", iEvt)); - - // Q2: how to get the SimCalorimeterHitCollection from the frame? (hint: check the collection name in the file) - const auto* simHits = static_cast(aframe.get("/*fill me*/")); - - float totalEnergy = 0.f; - - // loop over hits in the event - for (const auto& hit : *simHits) { - // Q3: how to decode the cellID to get the calo layer, abs layer, and sub-cell id? (hint: use the bit field - // decoder) - int caloLayer; /*fill me*/ - int absLayer; /*fill me*/ - int subCellId; /*fill me*/ - - // std::cout << "Hit energy: " << hit.getEnergy() << " GeV, CellID: " << hit.getCellID() - // << ", CaloLayer: " << caloLayer << ", AbsLayer: " << absLayer - // << ", SubCellId: " << subCellId << std::endl; - - // Q4: fill the histograms defined above to analyze the energy distribution, layer-wise energy sum - - // Q5: fill lateral shower shape per layer - // Note: cellid = 10*x + y - - // Q6: how to access the contributions to each hit and fill the timing distribution of contributions? (hint: check - // the edm4hep SimCalorimeterHit class definition for contributions) fill timing distribution of contributions - } // loop hits - } // loop events - - hLayerEnergySum->Scale(1. / static_cast(entries)); // average energy per layer - std::for_each(hLateralShapePerLayer.begin(), hLateralShapePerLayer.end(), [entries](TH2F* hist) { - hist->Scale(1. / static_cast(entries)); - }); // average lateral shape per layer - - // you know how to save or display the histograms from here - // output results - auto* outFile = TFile::Open("edm4hep_analysis.root", "RECREATE"); - hEnergy->Write(); - hLayerEnergySum->Write(); - hContributionTime->Write(); - - for (auto* hist : hLateralShapePerLayer) { - hist->Write(); - } - - outFile->Close(); -} - -// solution -// Q1: bit field decoder -// auto decoder = dd4hep::DDSegmentation::BitFieldCoder("calolayer:5,abslayer:1,cellid:10"); - -// Q2: get SimCalorimeterHitCollection from the frame -// const auto* simHits = static_cast(aframe.get("simplecaloRO")); - -// Q3: decode cellID -// int caloLayer = decoder.get(hit.getCellID(), "calolayer"); -// int absLayer = decoder.get(hit.getCellID(), "abslayer"); -// int subCellId = decoder.get(hit.getCellID(), "cellid"); - -// Q4: fill histograms -// total energy sum -// totalEnergy += hit.getEnergy(); -// // and add to energy histogram at the end of the event loop -// hEnergy->Fill(totalEnergy); - -// // energy per layer -// hLayerEnergySum->Fill(static_cast(caloLayer) + 0.5, hit.getEnergy()); -// // Note: the +0.5 is to fill the histogram bin corresponding to the integer layer number - -// Q5: fill lateral shower shape per layer -// int x = subCellId / 10; -// int y = subCellId % 10; -// unsigned int layerIndex = caloLayer - 1; // assuming caloLayer starts from 1 -// hLateralShapePerLayer[layerIndex]->Fill(static_cast(x) + 0.5, static_cast(y) + 0.5, hit.getEnergy()); - -// Q6: access contributions and fill timing distribution -// for (const auto& contrib : hit.getContributions()) { -// // std::cout << " Contrib energy: " << contrib.getEnergy() << " GeV, PDG: " << contrib.getPDG() << std::endl; -// hContributionTime->Fill(contrib.getTime(), contrib.getEnergy()); -// } diff --git a/DD4hepTutorials/sdaction/src/simplecaloSDAction.cpp b/DD4hepTutorials/sdaction/src/simplecaloSDAction.cpp index 09e33f0..db5af13 100644 --- a/DD4hepTutorials/sdaction/src/simplecaloSDAction.cpp +++ b/DD4hepTutorials/sdaction/src/simplecaloSDAction.cpp @@ -23,6 +23,7 @@ #include "DDG4/Geant4Mapping.h" #include "DDG4/Geant4SensDetAction.inl" +#include "CLHEP/Units/SystemOfUnits.h" #include "G4ThreeVector.hh" #include "G4TouchableHandle.hh" #include @@ -86,9 +87,12 @@ namespace sim { << "Vol " << aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetName() << " " << std::endl; #endif - dd4hep::BitFieldCoder decoder("calolayer:5,abslayer:1,cellid:10"); auto VolID = volumeID(aStep); #ifdef DEBUG + // Parsing the encoding string is not free, so keep it out of the per-step + // path: the decoder is only needed for the printout below. The string has to + // match the element of the block in simplecalo2.xml. + static const dd4hep::BitFieldCoder decoder("calolayer:5,abslayer:1,cellid:10"); auto CaloLayerID = decoder.get(VolID, "calolayer"); auto AbsLayerID = decoder.get(VolID, "abslayer"); auto CellID = decoder.get(VolID, "cellid"); @@ -98,35 +102,62 @@ namespace sim { G4TouchableHandle theTouchable = aStep->GetPreStepPoint()->GetTouchableHandle(); G4ThreeVector origin(0., 0., 0.); + // Centre of the cell this step happened in, in global coordinates G4ThreeVector CellPos = theTouchable->GetHistory()->GetTopTransform().Inverse().TransformPoint(origin); + // Where the step actually started, in global coordinates + G4ThreeVector StepPos = aStep->GetPreStepPoint()->GetPosition(); #ifdef DEBUG std::cout << "--> Cell global pos(mm) " << CellPos.x() << " " << CellPos.y() << " " << CellPos.z() << std::endl; #endif // Hands-on 5: apply a very short time cut (10 ns) to record the signals // and consider the cell border (2 cm) along x and y completely inefficient, - // i.e. not signal is recorder from that area. + // i.e. no signal is recorded from that area. // Hint: the x,y,z position of the step in the local volume reference frame is - // G4ThreeVector globalPosition = aStep->GetPreStepPoint()->GetPosition(); - // theTouchable->GetHistory()->GetTopTransform().TransformPoint(globalPosition); + // G4ThreeVector localPosition = + // theTouchable->GetHistory()->GetTopTransform().TransformPoint(StepPos); // // Hands-on 5: solution // /* - if (aStep->GetPreStepPoint()->GetGlobalTime() > 10) { + // Geant4 works in mm and ns, so CLHEP::ns and CLHEP::cm are the units to use + // here. Careful: inside namespace dd4hep, a bare `ns` or `cm` would resolve to + // dd4hep's own units, where the native length unit is cm and not mm. + constexpr double TIME_CUT = 10. * CLHEP::ns; + constexpr double CELL_HALF_SIZE = 5. * CLHEP::cm; // cells are 10 cm x 10 cm + constexpr double DEAD_BORDER = 2. * CLHEP::cm; // inefficient rim around each cell + constexpr double ACTIVE_HALF_SIZE = CELL_HALF_SIZE - DEAD_BORDER; + + if (aStep->GetPreStepPoint()->GetGlobalTime() > TIME_CUT) { return true; } - G4ThreeVector globalPosition = aStep->GetPreStepPoint()->GetPosition(); - G4ThreeVector localPosition = - theTouchable->GetHistory()->GetTopTransform().TransformPoint( - globalPosition); - if (std::abs(localPosition.x()) > 30. || std::abs(localPosition.y()) > 30.) { + G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(StepPos); + if (std::abs(localPosition.x()) > ACTIVE_HALF_SIZE || std::abs(localPosition.y()) > ACTIVE_HALF_SIZE) { return true; } // end of Hands-on 5 */ + // Describe this step as a contribution to the cell it happened in. + // + // Units: everything below stays in Geant4's internal units, i.e. MeV, mm and + // ns. Nothing is converted here, and nothing should be: DDG4's EDM4hep writer + // (Geant4Output2EDM4hep) divides by CLHEP::GeV, CLHEP::mm and CLHEP::ns when + // it serialises the hits, so the file ends up with the GeV, mm and ns that + // EDM4hep specifies. Converting here as well would apply the factor twice. + // Do not confuse either with dd4hep's own unit system, whose native length + // unit is cm. + Geant4Calorimeter::Hit::Contribution contrib; + contrib.trackID = aStep->GetTrack()->GetTrackID(); + contrib.pdgID = aStep->GetTrack()->GetParticleDefinition()->GetPDGEncoding(); + contrib.deposit = aStep->GetTotalEnergyDeposit(); + contrib.time = aStep->GetPreStepPoint()->GetGlobalTime(); + contrib.length = aStep->GetStepLength(); + // A contribution is a snapshot of a single step, so it carries the position of + // that step. EDM4hep stores this member as CaloHitContribution::stepPosition. + contrib.setPosition(StepPos.x(), StepPos.y(), StepPos.z()); + // Create the hits and accumulate contributions from multiple steps // Geant4HitCollection* coll = collection(m_collectionID); @@ -134,44 +165,18 @@ namespace sim { if (!hit) { // if the hit does not exist yet, create it hit = new Geant4Calorimeter::Hit(); - hit->cellID = VolID; // this should be assigned only once - // we divide the coordinated by 10 to save them as cm - Position HitCellPos(CellPos.x() / 10, CellPos.y() / 10, CellPos.z() / 10); - hit->position = HitCellPos; // this should be assigned only once - hit->energyDeposit = aStep->GetTotalEnergyDeposit(); - - // Add calo hit contributions - // - // Crete the first contribution associated to this hit - Geant4Calorimeter::Hit::Contribution contrib; - contrib.trackID = aStep->GetTrack()->GetTrackID(); - contrib.pdgID = aStep->GetTrack()->GetParticleDefinition()->GetPDGEncoding(); - contrib.deposit = aStep->GetTotalEnergyDeposit(); - contrib.time = aStep->GetPreStepPoint()->GetGlobalTime(); - contrib.x = HitCellPos.x(); - contrib.y = HitCellPos.y(); - contrib.z = HitCellPos.z(); - hit->truth.emplace_back(contrib); - + // cellID and position identify the readout element, so they are assigned + // exactly once: the hit sits at the centre of its cell. + hit->cellID = VolID; + hit->position = Position(CellPos.x(), CellPos.y(), CellPos.z()); + hit->energyDeposit = 0.; coll->add(VolID, hit); // add the hit to the hit collection - } else { // if the hit exists already, increment its fields - hit->energyDeposit += aStep->GetTotalEnergyDeposit(); - - // Add calo hit contributions - // - // Add a new contribution associated to this hit - Geant4Calorimeter::Hit::Contribution contrib; - contrib.trackID = aStep->GetTrack()->GetTrackID(); - contrib.pdgID = aStep->GetTrack()->GetParticleDefinition()->GetPDGEncoding(); - contrib.deposit = aStep->GetTotalEnergyDeposit(); - contrib.time = aStep->GetPreStepPoint()->GetGlobalTime(); - Position HitCellPos(CellPos.x() / 10, CellPos.y() / 10, CellPos.z() / 10); - contrib.x = HitCellPos.x(); - contrib.y = HitCellPos.y(); - contrib.z = HitCellPos.z(); - hit->truth.emplace_back(contrib); } + // Every step adds to the hit energy and appends its own contribution + hit->energyDeposit += aStep->GetTotalEnergyDeposit(); + hit->truth.emplace_back(contrib); + return true; } // end of Geant4SensitiveAction::process() method specialization diff --git a/DD4hepTutorials/simplecalo1/compact/simplecalo1.xml b/DD4hepTutorials/simplecalo1/compact/simplecalo1.xml index 6969c6e..f506dd3 100644 --- a/DD4hepTutorials/simplecalo1/compact/simplecalo1.xml +++ b/DD4hepTutorials/simplecalo1/compact/simplecalo1.xml @@ -95,6 +95,9 @@ grid_size_x="10*cm" grid_size_y="10*cm"/> + calolayer:5,abslayer:1,x:-10,y:-10 diff --git a/DD4hepTutorials/simplecalo1/sc1SteeringFile.py b/DD4hepTutorials/simplecalo1/sc1SteeringFile.py index 7dbf6e3..e7900ee 100644 --- a/DD4hepTutorials/simplecalo1/sc1SteeringFile.py +++ b/DD4hepTutorials/simplecalo1/sc1SteeringFile.py @@ -16,13 +16,22 @@ # See the License for the specific language governing permissions and # limitations under the License. # +import inspect +import os + +# Directory of this steering file, and the DD4hepTutorials directory above it, so +# that ddsim can be launched from anywhere. ddsim runs this file through exec() +# without defining __file__, hence inspect. +_HERE = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) +_TUTORIAL_DIR = os.path.dirname(_HERE) + from DDSim.DD4hepSimulation import DD4hepSimulation from g4units import mm, m, cm, GeV, MeV SIM = DD4hepSimulation() ## The compact XML file, or multiple compact files, if the last one is the closer. -SIM.compactFile = ["simplecalo1/compact/simplecalo1.xml"] +SIM.compactFile = [os.path.join(_HERE, "compact", "simplecalo1.xml")] ## Lorentz boost for the crossing angle, in radian! SIM.crossingAngleBoost = 0.0 SIM.enableDetailedShowerMode = False @@ -37,7 +46,10 @@ ## number of events to simulate, used in batch mode SIM.numberOfEvents = 500 ## Outputfile from the simulation: .slcio, edm4hep.root and .root output files are supported -SIM.outputFile = "simplecalo.root" +# Anchored on DD4hepTutorials/ rather than the current directory, because that is +# where the notebooks look for it (see drdcalo_tutorials.simplecalo1_input). +# --outputFile on the ddsim command line still overrides this. +SIM.outputFile = os.path.join(_TUTORIAL_DIR, "simplecalo1.root") ## Physics list to use in simulation SIM.physicsList = None ## Verbosity use integers from 1(most) to 7(least) verbose diff --git a/DD4hepTutorials/simplecalo2/include/sc2_solution1.h b/DD4hepTutorials/simplecalo2/include/sc2_solution1.h index 9601011..1685574 100644 --- a/DD4hepTutorials/simplecalo2/include/sc2_solution1.h +++ b/DD4hepTutorials/simplecalo2/include/sc2_solution1.h @@ -22,12 +22,30 @@ CellVol.setVisAttributes(description, x_cell.visStr()); // Make the cell sensitive if (iscellsens) CellVol.setSensitiveDetector(sens); -double x, y = 0; -for (std::size_t i = 0; i < 10; i++) { - y = SensLayerY / 2. - CellY / 2. - i * CellY; - for (std::size_t j = 0; j < 10; j++) { - x = -SensLayerX / 2. + CellX / 2. + j * CellX; - PlacedVolume CellVolPlaced = SensLayerVol.placeVolume(CellVol, 10 * i + j, Position(x, y, 0.)); - CellVolPlaced.addPhysVolID("cellid", 10 * i + j); + +// How many cells fit into the sensitive layer follows from the XML dimensions, so +// changing CellX/CellY or SensLayerX/SensLayerY in simplecalo2.xml is enough and +// this code does not have to be touched. +const int NCellsX = static_cast(std::round(SensLayerX / CellX)); +const int NCellsY = static_cast(std::round(SensLayerY / CellY)); + +// Refuse a geometry that would be silently wrong. Cells that do not tile the layer +// leave part of it uncovered, and the cellid field of the readout is 10 bits wide +// (see in simplecalo2.xml), so beyond 1024 cells two of them share an address. +if (std::abs(NCellsX * CellX - SensLayerX) > 1e-6 || std::abs(NCellsY * CellY - SensLayerY) > 1e-6) + except("simplecalo2", "Cells do not tile the sensitive layer: %g/%g and %g/%g must be whole numbers.", SensLayerX, + CellX, SensLayerY, CellY); +if (NCellsX * NCellsY > 1024) + except("simplecalo2", "%d x %d cells do not fit the 10-bit cellid field, which holds 1024.", NCellsX, NCellsY); + +// The outer loop runs over y and the inner one over x, hence cellid = NCellsX * iY + iX. +// The analysis has to use the same convention to turn a cellid back into a cell position. +for (int iY = 0; iY < NCellsY; iY++) { + const double y = SensLayerY / 2. - CellY / 2. - iY * CellY; + for (int iX = 0; iX < NCellsX; iX++) { + const double x = -SensLayerX / 2. + CellX / 2. + iX * CellX; + const int cellid = NCellsX * iY + iX; + PlacedVolume CellVolPlaced = SensLayerVol.placeVolume(CellVol, cellid, Position(x, y, 0.)); + CellVolPlaced.addPhysVolID("cellid", cellid); } } diff --git a/DD4hepTutorials/simplecalo2/sc2SteeringFile.py b/DD4hepTutorials/simplecalo2/sc2SteeringFile.py index 86ce0c8..39fed41 100644 --- a/DD4hepTutorials/simplecalo2/sc2SteeringFile.py +++ b/DD4hepTutorials/simplecalo2/sc2SteeringFile.py @@ -16,13 +16,22 @@ # See the License for the specific language governing permissions and # limitations under the License. # +import inspect +import os + +# Directory of this steering file, and the DD4hepTutorials directory above it, so +# that ddsim can be launched from anywhere. ddsim runs this file through exec() +# without defining __file__, hence inspect. +_HERE = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) +_TUTORIAL_DIR = os.path.dirname(_HERE) + from DDSim.DD4hepSimulation import DD4hepSimulation from g4units import mm, m, cm, GeV, MeV SIM = DD4hepSimulation() ## The compact XML file, or multiple compact files, if the last one is the closer. -SIM.compactFile = ["simplecalo2/compact/simplecalo2.xml"] +SIM.compactFile = [os.path.join(_HERE, "compact", "simplecalo2.xml")] ## Lorentz boost for the crossing angle, in radian! SIM.crossingAngleBoost = 0.0 SIM.enableDetailedShowerMode = False @@ -37,7 +46,10 @@ ## number of events to simulate, used in batch mode SIM.numberOfEvents = 500 ## Outputfile from the simulation: .slcio, edm4hep.root and .root output files are supported -SIM.outputFile = "simplecalo.root" +# Anchored on DD4hepTutorials/ rather than the current directory, because that is +# where the notebooks look for it (see drdcalo_tutorials.simplecalo2_input). +# --outputFile on the ddsim command line still overrides this. +SIM.outputFile = os.path.join(_TUTORIAL_DIR, "simplecalo2.root") ## Physics list to use in simulation SIM.physicsList = None ## Verbosity use integers from 1(most) to 7(least) verbose diff --git a/DD4hepTutorials/simplecalo2/src/simplecalo2.cpp b/DD4hepTutorials/simplecalo2/src/simplecalo2.cpp index 8f7428c..b811edb 100644 --- a/DD4hepTutorials/simplecalo2/src/simplecalo2.cpp +++ b/DD4hepTutorials/simplecalo2/src/simplecalo2.cpp @@ -20,6 +20,8 @@ // Includers from DD4hep #include "DDRec/Vector3D.h" #include +#include // dd4hep::except, used by sc2_solution1.h +#include // std::round and std::abs, likewise using namespace dd4hep; diff --git a/GaudiTutorial/EventStats/options/runEventStatsMultithreading.py b/GaudiTutorial/EventStats/options/runEventStatsMultithreading.py index a168fe7..4da9955 100644 --- a/GaudiTutorial/EventStats/options/runEventStatsMultithreading.py +++ b/GaudiTutorial/EventStats/options/runEventStatsMultithreading.py @@ -16,6 +16,8 @@ # See the License for the specific language governing permissions and # limitations under the License. # +from drdcalo_tutorials import gaudi_output, sample_input + from Gaudi.Configuration import INFO from k4FWCore import IOSvc, ApplicationMgr from Configurables import EventDataSvc, AuditorSvc, ChronoAuditor, HiveWhiteBoard, HiveSlimEventLoopMgr, AvalancheSchedulerSvc @@ -27,8 +29,10 @@ scheduler = AvalancheSchedulerSvc(ThreadPoolSize=threads) io_svc = IOSvc("IOSvc") -io_svc.Input = "../../data/simpleCalo_simulation.root" -io_svc.Output = "../../data/simpleCalo_eventStats.root" +io_svc.Input = sample_input() +# Its own output names, so that running this next to runEventStatsSolution.py +# compares the two results instead of overwriting one with the other. +io_svc.Output = gaudi_output("simpleCalo_eventStats_multithreaded.root") chra = ChronoAuditor() audsvc = AuditorSvc() @@ -41,6 +45,7 @@ OutputEnergyBarycentre = ["EnergyBarycentreX", "EnergyBarycentreY", "EnergyBarycentreZ"], OutputTotalEnergy = ["TotalEnergy"], SaveHistograms = True, + HistogramFile = gaudi_output("debugEventStats_multithreaded.root"), OutputLevel = INFO ) diff --git a/GaudiTutorial/EventStats/options/runEventStatsSolution.py b/GaudiTutorial/EventStats/options/runEventStatsSolution.py index f0b425d..a7bcee0 100644 --- a/GaudiTutorial/EventStats/options/runEventStatsSolution.py +++ b/GaudiTutorial/EventStats/options/runEventStatsSolution.py @@ -16,13 +16,15 @@ # See the License for the specific language governing permissions and # limitations under the License. # +from drdcalo_tutorials import gaudi_output, sample_input + from Gaudi.Configuration import INFO from k4FWCore import IOSvc, ApplicationMgr from Configurables import EventDataSvc, AuditorSvc, ChronoAuditor io_svc = IOSvc("IOSvc") -io_svc.Input = "../../data/simpleCalo_simulation.root" -io_svc.Output = "../../data/simpleCalo_eventStats.root" +io_svc.Input = sample_input() +io_svc.Output = gaudi_output("simpleCalo_eventStats.root") chra = ChronoAuditor() audsvc = AuditorSvc() @@ -35,6 +37,7 @@ OutputEnergyBarycentre = ["EnergyBarycentreX", "EnergyBarycentreY", "EnergyBarycentreZ"], OutputTotalEnergy = ["TotalEnergy"], SaveHistograms = True, + HistogramFile = gaudi_output("debugEventStats.root"), OutputLevel = INFO ) diff --git a/GaudiTutorial/EventStats/src/components/EventStats.cpp b/GaudiTutorial/EventStats/src/components/EventStats.cpp index 15e4d94..2005728 100644 --- a/GaudiTutorial/EventStats/src/components/EventStats.cpp +++ b/GaudiTutorial/EventStats/src/components/EventStats.cpp @@ -27,6 +27,8 @@ // STL #include +#include +#include #include #include @@ -134,7 +136,6 @@ struct EventStats final barycentre_x += x * energy; barycentre_y += y * energy; - barycentre_y += y * energy; barycentre_z += z * energy; totalEnergy += energy; @@ -164,15 +165,23 @@ struct EventStats final min_pos_z = z; } - barycentre_x /= totalEnergy; - barycentre_y /= totalEnergy; - barycentre_z /= totalEnergy; + // An event without any energy deposit has no barycentre: dividing here would + // hand NaN to every downstream algorithm. + if (totalEnergy > 0.0) { + barycentre_x /= totalEnergy; + barycentre_y /= totalEnergy; + barycentre_z /= totalEnergy; + } ///////////////////////////////////////// //// SAVE HISTOGRAMS OF ENERGY STATS //// ///////////////////////////////////////// - if (m_saveHisto) { + if (m_saveHisto && totalEnergy > 0.0) { + // operator() is const and the Avalanche scheduler calls it from several + // threads at once (see runEventStatsMultithreading.py). TH1::Fill is not + // thread safe, so the shared histograms have to be protected. + std::scoped_lock lock(m_histoMutex); hTotalEnergy->Fill(totalEnergy); hMaxEnergy->Fill(maxEnergy); hMinEnergy->Fill(minEnergy); @@ -199,14 +208,28 @@ struct EventStats final info() << "Event statistics:" << endmsg; info() << "Total deposited energy = " << totalEnergy << " GeV" << endmsg; - info() << "Energy Barycentre position: (" << barycentre_x << ", " << barycentre_y << ", " << barycentre_z << ") mm" - << endmsg; - info() << "Max energy: " << maxEnergy << " GeV at (" << maxEnergy_x << ", " << maxEnergy_y << ", " << maxEnergy_z - << ") mm" << endmsg; - info() << "Min energy: " << minEnergy << " GeV at (" << minEnergy_x << ", " << minEnergy_y << ", " << minEnergy_z - << ") mm" << endmsg; - info() << "MaxXYZ values: (" << max_pos_x << ", " << max_pos_y << ", " << max_pos_z << ") mm" << endmsg; - info() << "MinXYZ values: (" << min_pos_x << ", " << min_pos_y << ", " << min_pos_z << ") mm" << endmsg; + + // An event whose collection is empty is a normal occurrence: the particle can + // miss the calorimeter. The loop above then never ran, so every extremum is + // still at the sentinel it was initialised with, and printing + // numeric_limits::lowest() as "the maximum energy" is worse than saying + // nothing at all. + if (InputSimCaloHitCollection.size() == 0) { + info() << "No calorimeter hits in this event: no barycentre and no extrema." << endmsg; + } else { + if (totalEnergy > 0.0) { + info() << "Energy Barycentre position: (" << barycentre_x << ", " << barycentre_y << ", " << barycentre_z + << ") mm" << endmsg; + } else { + info() << "No energy deposited in this event: the barycentre is undefined." << endmsg; + } + info() << "Max energy: " << maxEnergy << " GeV at (" << maxEnergy_x << ", " << maxEnergy_y << ", " << maxEnergy_z + << ") mm" << endmsg; + info() << "Min energy: " << minEnergy << " GeV at (" << minEnergy_x << ", " << minEnergy_y << ", " << minEnergy_z + << ") mm" << endmsg; + info() << "MaxXYZ values: (" << max_pos_x << ", " << max_pos_y << ", " << max_pos_z << ") mm" << endmsg; + info() << "MinXYZ values: (" << min_pos_x << ", " << min_pos_y << ", " << min_pos_z << ") mm" << endmsg; + } info() << "------------------------------" << endmsg; //////////////// @@ -221,7 +244,7 @@ struct EventStats final if (m_saveHisto.value()) { - TFile* f = new TFile("debugEventStats.root", "RECREATE"); + TFile* f = new TFile(m_histoFile.value().c_str(), "RECREATE"); hTotalEnergy->Write(); hMaxEnergy->Write(); hMinEnergy->Write(); @@ -239,9 +262,16 @@ struct EventStats final private: Gaudi::Property m_saveHisto{this, "SaveHistograms", false, "flag to save histograms"}; - TH1D* hTotalEnergy; - TH1D* hMaxEnergy; - TH1D* hMinEnergy; + // Written relative to the current directory, so the tests point it at their own + // build directory rather than dropping a file into the source tree. + Gaudi::Property m_histoFile{this, "HistogramFile", "debugEventStats.root", + "file the histograms above are written to"}; + + mutable std::mutex m_histoMutex; + + TH1D* hTotalEnergy = nullptr; + TH1D* hMaxEnergy = nullptr; + TH1D* hMinEnergy = nullptr; }; DECLARE_COMPONENT(EventStats) diff --git a/GaudiTutorial/MoliereRadiusFunctional/options/runMoliereRadius.py b/GaudiTutorial/MoliereRadiusFunctional/options/runMoliereRadius.py index 11d0b5f..2a09577 100644 --- a/GaudiTutorial/MoliereRadiusFunctional/options/runMoliereRadius.py +++ b/GaudiTutorial/MoliereRadiusFunctional/options/runMoliereRadius.py @@ -16,6 +16,8 @@ # See the License for the specific language governing permissions and # limitations under the License. # +from drdcalo_tutorials import gaudi_output, sample_input + from Gaudi.Configuration import INFO from k4FWCore import IOSvc, ApplicationMgr from Configurables import EventDataSvc, UniqueIDGenSvc, ChronoAuditor, AuditorSvc @@ -25,8 +27,8 @@ audsvc.Auditors = [chra] io_svc = IOSvc("IOSvc") -io_svc.Input = "../../data/simpleCalo_simulation.root" -io_svc.Output = "../../data/simpleCalo_moliereRadius.root" +io_svc.Input = sample_input() +io_svc.Output = gaudi_output("simpleCalo_moliereRadius.root") from Configurables import EventStats diff --git a/GaudiTutorial/MoliereRadiusFunctional/options/runMoliereRadiusSolution.py b/GaudiTutorial/MoliereRadiusFunctional/options/runMoliereRadiusSolution.py index b94c66e..2f529ac 100644 --- a/GaudiTutorial/MoliereRadiusFunctional/options/runMoliereRadiusSolution.py +++ b/GaudiTutorial/MoliereRadiusFunctional/options/runMoliereRadiusSolution.py @@ -16,6 +16,8 @@ # See the License for the specific language governing permissions and # limitations under the License. # +from drdcalo_tutorials import gaudi_output, sample_input + from Gaudi.Configuration import INFO from k4FWCore import IOSvc, ApplicationMgr from Configurables import EventDataSvc, UniqueIDGenSvc, ChronoAuditor, AuditorSvc @@ -25,8 +27,8 @@ audsvc.Auditors = [chra] io_svc = IOSvc("IOSvc") -io_svc.Input = "../../data/simpleCalo_simulation.root" -io_svc.Output = "../../data/simpleCalo_moliereRadius.root" +io_svc.Input = sample_input() +io_svc.Output = gaudi_output("simpleCalo_moliereRadius.root") from Configurables import EventStats diff --git a/GaudiTutorial/MoliereRadiusFunctional/src/components/MoliereRadiusSolution.cpp b/GaudiTutorial/MoliereRadiusFunctional/src/components/MoliereRadiusSolution.cpp index 6a3d79e..e85a42d 100644 --- a/GaudiTutorial/MoliereRadiusFunctional/src/components/MoliereRadiusSolution.cpp +++ b/GaudiTutorial/MoliereRadiusFunctional/src/components/MoliereRadiusSolution.cpp @@ -30,6 +30,7 @@ #include "edm4hep/CalorimeterHitCollection.h" // STL +#include // std::sort #include #include #include diff --git a/GaudiTutorial/README.md b/GaudiTutorial/README.md index df230e7..5eba335 100644 --- a/GaudiTutorial/README.md +++ b/GaudiTutorial/README.md @@ -16,7 +16,7 @@ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. --> -# GaudiTutorials for DRD6 +# GaudiTutorials for DRDCalo This repository hosts tutorials for using the *Gaudi* software as part of the larger *key4hep* ecosystem. @@ -24,15 +24,21 @@ It contains hands-on exercises to get familiar with Gaudi steering files and alg To experiment with the tutorial, you can follow [this presentation](https://indico.cern.ch/event/1551941/sessions/613344/#20250918) at the 4th DRD-Calo Collaboration Meeting, which covers the *EventStats*, *RandomNoiseDigitizer*, and *MoliereRadius* exercises. This README contains a short description for each of the hosted exercises. -The exercises run on data that has been created with with the *simplecalo* calorimeter from the **DD4hepTutorials** exercises. -An example data file will be downloaded automatically when compiling the repository. +The exercises run on data that has been created with the *simplecalo* calorimeter from the **DD4hepTutorials** exercises. +No data download is required: the steering files call +`drdcalo_tutorials.sample_input()`, the bundled 10-event +`DD4hepTutorials/data/simplecalo2_sample.root`, and write through +`drdcalo_tutorials.gaudi_output(".root")` into `GaudiTutorial/data/`. Both can be +overridden on the command line with `--IOSvc.Input` and `--IOSvc.Output`, so a full +`simplecalo2.root` from the DD4hep tutorial can be analysed with +`k4run --IOSvc.Input DD4hepTutorials/simplecalo2.root`. ## EventStats The goal of this exercise is to become familiar with the Gaudi steering file. For this purpose, an EventStats algorithm is provided, which saves the energy barycentre and total energy for each event using the `podio::UserDataCollection`. -You should adapt the steering file `runEventStats.py` in the `EventStats/options' folder such that this algorithm is executed on the data. +You should adapt the steering file `runEventStats.py` in the `EventStats/options` folder such that this algorithm is executed on the data. A solution file is provided. ## RandomNoiseDigitizer diff --git a/GaudiTutorial/RandomNoiseDigitizer/options/runRandomNoiseDigitizer.py b/GaudiTutorial/RandomNoiseDigitizer/options/runRandomNoiseDigitizer.py index bbef440..c26a93c 100644 --- a/GaudiTutorial/RandomNoiseDigitizer/options/runRandomNoiseDigitizer.py +++ b/GaudiTutorial/RandomNoiseDigitizer/options/runRandomNoiseDigitizer.py @@ -16,6 +16,8 @@ # See the License for the specific language governing permissions and # limitations under the License. # +from drdcalo_tutorials import gaudi_output, sample_input + from Gaudi.Configuration import INFO from k4FWCore import IOSvc, ApplicationMgr from Configurables import EventDataSvc, UniqueIDGenSvc, ChronoAuditor, AuditorSvc @@ -25,8 +27,8 @@ audsvc.Auditors = [chra] io_svc = IOSvc("IOSvc") -io_svc.Input = "../../data/simpleCalo_simulation.root" -io_svc.Output = "../../data/simpleCalo_noiseDigitizer.root" +io_svc.Input = sample_input() +io_svc.Output = gaudi_output("simpleCalo_noiseDigitizer.root") from Configurables import EventStats diff --git a/GaudiTutorial/RandomNoiseDigitizer/options/runRandomNoiseDigitizerSolution.py b/GaudiTutorial/RandomNoiseDigitizer/options/runRandomNoiseDigitizerSolution.py index bc09985..cbff590 100644 --- a/GaudiTutorial/RandomNoiseDigitizer/options/runRandomNoiseDigitizerSolution.py +++ b/GaudiTutorial/RandomNoiseDigitizer/options/runRandomNoiseDigitizerSolution.py @@ -16,6 +16,8 @@ # See the License for the specific language governing permissions and # limitations under the License. # +from drdcalo_tutorials import gaudi_output, sample_input + from Gaudi.Configuration import INFO from k4FWCore import IOSvc, ApplicationMgr from Configurables import EventDataSvc, UniqueIDGenSvc, ChronoAuditor, AuditorSvc @@ -25,8 +27,8 @@ audsvc.Auditors = [chra] io_svc = IOSvc("IOSvc") -io_svc.Input = "../../data/simpleCalo_simulation.root" -io_svc.Output = "../../data/simpleCalo_noiseDigitizer.root" +io_svc.Input = sample_input() +io_svc.Output = gaudi_output("simpleCalo_noiseDigitizer.root") from Configurables import EventStats diff --git a/README.md b/README.md index efab6cb..3397c6a 100644 --- a/README.md +++ b/README.md @@ -1,32 +1,82 @@ # SoftwareTutorials - -This repository hosts a variety of software tutorials related to DRDCalo (formerly DRD6) Collaboriation activities. -It includes the tutorials on DD4hep and Gaudi, presented in the 2nd and 4th Collaboration Meeting respectively. - -Ideally, this repository is expanded with relevant tutorials in the future. - -The tutorials can be completed by following the presentation slides linked in the corresponding sub-directories - +This repository hosts software tutorials related to DRDCalo (formerly DRD6) Collaboration +activities: DD4hep and Gaudi, presented at the 2nd and 4th Collaboration Meeting respectively. +Each is completed by following the presentation slides linked in its sub-directory. ## Compilation -These tutorials are made within the key4hep environment. -To compile and run them, access to an Alma9 machine with `/cvmfs/` mounted (e.g. lxplus) to source the key4hep stack is required. - -To clone and build the repository, run the following commands: - +These tutorials run inside the key4hep environment, so they need an EL9 machine (AlmaLinux 9, +RHEL 9, ...) with `/cvmfs/` mounted — lxplus, for instance. ``` bash -source /cvmfs/sw.hsf.org/key4hep/setup.sh -git clone https://github.com/DRD6/SoftwareTutorials.git +git clone https://github.com/DRDCalo/SoftwareTutorials.git cd SoftwareTutorials -k4_local_repo +source setup.sh mkdir build install cd build cmake .. -DCMAKE_INSTALL_PREFIX=../install make install -j6 ``` -The `k4_local_repo` command sets the environment variables and makes sure the local modules are found by the system. -It needs to be run in **every new shell**. +`setup.sh` must be **sourced in every new shell**. It + +1. sources the key4hep stack, pinned to a fixed release so that everybody works with the same + software (override with `KEY4HEP_VERSION=... source setup.sh`); +2. runs `k4_local_repo`, pointing the environment at the libraries built here; +3. creates `.venv`, a virtual environment inheriting everything from the stack; +4. makes `drdcalo_tutorials` importable — the module resolving every tutorial's input and output + paths — through `PYTHONPATH` and a `.pth` file in `.venv`, so even a process that does not + inherit the variable finds it; +5. writes `.env` with the key4hep runtime paths, which is what an editor hands to a notebook + kernel. Only those variables are written, never the rest of your shell environment, and the + file is created `0600`. + +If a different key4hep release is already active, `setup.sh` stops rather than mixing the two. +Inside a container step 5 is skipped and it says so: `.env` records absolute paths and is read by +an editor on the host, where the container's mount point would be wrong. + +Steps 3 to 5 exist for the notebooks: an editor's Jupyter kernel does not inherit your shell, so +it needs an interpreter it can find (`.venv`) and the environment to run in (`.env`). Both are +ignored by git and are rebuilt by sourcing `setup.sh` again. + +### Using the notebooks in VS Code + +**Open `SoftwareTutorials` itself as the workspace folder**, not a parent directory. This is the +one thing that has to be right: VS Code applies `${workspaceFolder}/.env`, and from a level above +it finds nothing and the notebook fails with `ModuleNotFoundError: No module named 'awkward'`. +Then choose **Select Kernel** → **Python Environments** and pick `.venv/bin/python` (shown as +*DRDCalo Tutorial*). No kernel has to be registered. + +To keep a parent directory as the workspace folder instead, point VS Code at the right file in +that workspace's `.vscode/settings.json`: + +``` json +{ + "python.envFile": "${workspaceFolder}/SoftwareTutorials/.env" +} +``` + +If a module goes missing, the recorded environment is stale: source `setup.sh` again and restart +the kernel. If it stays missing, the variable providing it is not on the allow-list in step 5 of +`setup.sh` — add the name there. If VS Code offers to install packages, it is running an +interpreter without the key4hep stack; fix the workspace folder rather than accepting. + +### Disk space + +`ddsim` writes a large EDM4hep file: the 500 events configured in the steering files are roughly +600 MB, and each tutorial writes its own. Without `--outputFile` they land in +`DD4hepTutorials/simplecalo1.root` and `DD4hepTutorials/simplecalo2.root` — where the notebooks +look — whatever directory `ddsim` was launched from. On a small quota, write elsewhere and point +the notebook at it: + +``` bash +ddsim --steeringFile DD4hepTutorials/simplecalo1/sc1SteeringFile.py --outputFile /tmp/simplecalo1.root +export SIMPLECALO1_FILE=/tmp/simplecalo1.root +``` + +or lower `SIM.numberOfEvents`; 100 events are already plenty for the energy resolution fit. + +Nothing forces you to generate that file at all: the notebooks and all Gaudi exercises fall back +to the bundled 10-event `DD4hepTutorials/data/simplecalo2_sample.root`, and say so when they do. +Every one of these rules lives in `python/drdcalo_tutorials/__init__.py`. diff --git a/python/drdcalo_tutorials/__init__.py b/python/drdcalo_tutorials/__init__.py new file mode 100644 index 0000000..927ea8b --- /dev/null +++ b/python/drdcalo_tutorials/__init__.py @@ -0,0 +1,142 @@ +# +# Copyright (c) 2020-2024 Key4hep-Project. +# +# This file is part of Key4hep. +# See https://key4hep.github.io/key4hep-doc/ for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +"""Where the tutorials read their input and write their output. + +Every exercise resolves its files through this module, so the layout of the +repository is written down once instead of once per script. `setup.sh` puts +this directory on PYTHONPATH, which is also how the tutorial plugins become +importable, so anything that can run an exercise can import this. + +The module sits at /python/drdcalo_tutorials, which is what makes +`REPOSITORY` below correct. The helpers return strings, ready to be assigned to +a Gaudi property. +""" + +import os +import sys +from pathlib import Path +from xml.etree import ElementTree + +REPOSITORY = Path(__file__).resolve().parents[2] + +DD4HEP_TUTORIALS = REPOSITORY / "DD4hepTutorials" +GAUDI_TUTORIALS = REPOSITORY / "GaudiTutorial" + +#: Generated Gaudi output. Kept out of the source directories, and out of git. +GAUDI_DATA = GAUDI_TUTORIALS / "data" + +#: The 10-event sample committed to the repository, so that the exercises run +#: straight after a clone without a multi-hundred-MB simulation. +SAMPLE_INPUT = DD4HEP_TUTORIALS / "data" / "simplecalo2_sample.root" + +#: Where the simplecalo2 steering file writes its full simulation output. +SIMULATED_INPUT = DD4HEP_TUTORIALS / "simplecalo2.root" + +#: Where the simplecalo1 steering file writes its full simulation output. +SIMULATED_INPUT_SC1 = DD4HEP_TUTORIALS / "simplecalo1.root" + +#: The compact files the two exercises build their geometry from. Hands-on 6 reads +#: its cell dimensions back out of the simplecalo2 one through compact_constants() +#: instead of repeating the numbers, so that the analysis follows the geometry +#: whenever the XML is changed. +SIMPLECALO1_COMPACT = DD4HEP_TUTORIALS / "simplecalo1" / "compact" / "simplecalo1.xml" +SIMPLECALO2_COMPACT = DD4HEP_TUTORIALS / "simplecalo2" / "compact" / "simplecalo2.xml" + + +def sample_input() -> str: + """The bundled 10-event sample. + + Used by the Gaudi exercises: it is small, quick, and there right after a + clone, so nobody has to run the simulation first. Pass --IOSvc.Input to + k4run to analyse a different file. + """ + return str(SAMPLE_INPUT) + + +def _announce_fallback(expected: Path, steering: str) -> None: + """Say out loud that the bundled 10-event sample is being used. + + Silently analysing a different file than the one the reader thinks they + produced is the most confusing thing this module could do, so the fallback + is never quiet. stderr, because in a notebook that is the stream that stands + out from the analysis output. + """ + print( + f"NOTE: {expected} does not exist, falling back to the bundled 10-event sample.\n" + f" Ten events are enough to make the code run, not to make a plot worth " + f"showing.\n" + f" Produce the full simulation with: ddsim --steeringFile {steering}", + file=sys.stderr, + ) + + +def simplecalo2_input() -> str: + """The simplecalo2 events to analyse in the notebooks. + + A full simulation generated in the DD4hep tutorial is preferred when it is + there; otherwise the bundled sample is used. Set SIMPLECALO2_FILE to read + another compatible file instead. Note that the full simulation is several + hundred MB, which is why the Gaudi exercises stay on sample_input(). + """ + override = os.getenv("SIMPLECALO2_FILE") + if override: + return str(Path(override).expanduser()) + if SIMULATED_INPUT.exists(): + return str(SIMULATED_INPUT) + _announce_fallback(SIMULATED_INPUT, "simplecalo2/sc2SteeringFile.py") + return str(SAMPLE_INPUT) + + +def simplecalo1_input() -> str: + """The simplecalo1 events to analyse in the Section 1 notebook. + + A full simulation generated in the DD4hep tutorial is preferred when it is + there; otherwise the bundled sample is used, which carries the same + `simplecaloRO` collection and so gives a meaningful energy sum even though it + was produced with simplecalo2. Set SIMPLECALO1_FILE to read another file. + """ + override = os.getenv("SIMPLECALO1_FILE") + if override: + return str(Path(override).expanduser()) + if SIMULATED_INPUT_SC1.exists(): + return str(SIMULATED_INPUT_SC1) + _announce_fallback(SIMULATED_INPUT_SC1, "simplecalo1/sc1SteeringFile.py") + return str(SAMPLE_INPUT) + + +def gaudi_output(filename: str) -> str: + """Path for a file written by a Gaudi exercise, creating the directory.""" + GAUDI_DATA.mkdir(parents=True, exist_ok=True) + return str(GAUDI_DATA / filename) + + +def compact_constants(compact) -> dict: + """The definitions of a compact file, lengths in mm. + + `` comes back as `{"CellX": 100.0}`. + Values are evaluated in the order the file defines them, so a constant can + refer to an earlier one the way the compact file does. mm is also the unit + EDM4hep stores positions in, so what comes back compares with a hit position + directly. + """ + values = {"mm": 1.0, "cm": 10.0, "m": 1000.0} + for constant in ElementTree.parse(str(compact)).iter("constant"): + name, value = constant.get("name"), constant.get("value") + values[name] = eval(value, {"__builtins__": {}}, values) + return values diff --git a/setup.sh b/setup.sh new file mode 100755 index 0000000..b79807f --- /dev/null +++ b/setup.sh @@ -0,0 +1,279 @@ +#!/bin/bash +## +## Copyright (c) 2020-2024 Key4hep-Project. +## +## This file is part of Key4hep. +## See https://key4hep.github.io/key4hep-doc/ for further info. +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. +## + +# This script has to be sourced: executing it would set up a shell that exits +# immediately afterwards. +if ! (return 0 2>/dev/null); then + echo "This script must be sourced: source setup.sh" + exit 1 +fi + +_TUTORIAL_ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +_TUTORIAL_PREV_DIR="${PWD}" + +# 1. key4hep stack, pinned so that everybody follows the tutorial on the same +# software. Reuse it only when the requested release is already active. +KEY4HEP_VERSION="${KEY4HEP_VERSION:-2026-04-08}" +_TUTORIAL_NEED_STACK=true +case "${KEY4HEP_STACK:-}" in + */releases/"${KEY4HEP_VERSION}"/*) + if command -v k4run > /dev/null 2>&1; then + echo "key4hep ${KEY4HEP_VERSION} already set up, skipping" + _TUTORIAL_NEED_STACK=false + else + # KEY4HEP_STACK is exported but the paths it describes are not in this + # shell. That happens whenever a process inherits the variable without + # the environment, for instance a terminal opened by a VS Code server + # that was itself started from a key4hep shell. The stack refuses to set + # itself up again while the variable is there, so clear it first. + echo "key4hep ${KEY4HEP_VERSION} is recorded but not active here, setting it up again ..." + unset KEY4HEP_STACK + fi + ;; + "") + echo "Setting up key4hep ${KEY4HEP_VERSION} ..." + ;; + *) + echo "A different key4hep release is already active: ${KEY4HEP_STACK}" + echo "Start a fresh shell without it, then source setup.sh again." + return 1 + ;; +esac + +if ${_TUTORIAL_NEED_STACK}; then + if ! source /cvmfs/sw.hsf.org/key4hep/setup.sh -r "${KEY4HEP_VERSION}"; then + echo "Failed to set up key4hep ${KEY4HEP_VERSION}." + return 1 + fi +fi + +# The stack prints rather than fails in a few situations, so confirm that its +# tools are really reachable before building anything on top of it. +if ! command -v k4run > /dev/null 2>&1; then + echo "key4hep ${KEY4HEP_VERSION} did not set up correctly: k4run is not on PATH." + echo "Start a fresh shell and source setup.sh again." + return 1 +fi + +# 2. Make the libraries built in this repository visible. +cd "${_TUTORIAL_ROOT}" || return 1 +k4_local_repo +cd "${_TUTORIAL_PREV_DIR}" || return 1 + +# 3. A virtual environment that inherits everything from the stack. Repair an +# incomplete environment too; VS Code can otherwise leave behind an empty +# .venv after a failed attempt to install notebook dependencies. +if [ ! -x "${_TUTORIAL_ROOT}/.venv/bin/python3" ] || \ + [ ! -f "${_TUTORIAL_ROOT}/.venv/bin/activate" ] || \ + ! grep -q '^include-system-site-packages = true$' \ + "${_TUTORIAL_ROOT}/.venv/pyvenv.cfg" 2>/dev/null; then + if [ -d "${_TUTORIAL_ROOT}/.venv" ]; then + echo "Repairing incomplete .venv ..." + else + echo "Creating .venv ..." + fi + python -m venv --clear "${_TUTORIAL_ROOT}/.venv" \ + --system-site-packages --prompt "DRDCalo Tutorial" +fi + +source "${_TUTORIAL_ROOT}/.venv/bin/activate" + +# venv activation scripts contain the absolute path used when they were +# created. Normalise it here so that a moved checkout (or a checkout created +# through a mounted workspace) still selects this repository's interpreter. +export VIRTUAL_ENV="${_TUTORIAL_ROOT}/.venv" + +# Move the environment to the front unconditionally rather than only adding it +# when absent: step 1 prepends the stack's own bin directories, so on a second +# `source setup.sh` in the same shell the stack python would otherwise sit ahead +# of this one and the rest of this script would configure the wrong interpreter. +_TUTORIAL_PATH="" +_TUTORIAL_IFS="${IFS}" +IFS=':' +for _TUTORIAL_ENTRY in ${PATH}; do + [ "${_TUTORIAL_ENTRY}" = "${VIRTUAL_ENV}/bin" ] && continue + _TUTORIAL_PATH="${_TUTORIAL_PATH:+${_TUTORIAL_PATH}:}${_TUTORIAL_ENTRY}" +done +IFS="${_TUTORIAL_IFS}" +export PATH="${VIRTUAL_ENV}/bin${_TUTORIAL_PATH:+:${_TUTORIAL_PATH}}" +unset _TUTORIAL_PATH _TUTORIAL_IFS _TUTORIAL_ENTRY +hash -r 2>/dev/null + +# 4. Make the drdcalo_tutorials module importable. Every exercise resolves its +# input and output files through it, so it has to be on PYTHONPATH before +# the environment is captured below. +case ":${PYTHONPATH:-}:" in + *":${_TUTORIAL_ROOT}/python:"*) ;; + *) export PYTHONPATH="${_TUTORIAL_ROOT}/python${PYTHONPATH:+:${PYTHONPATH}}" ;; +esac + +# A .pth file in the virtual environment makes the module importable by +# .venv/bin/python whatever PYTHONPATH it is started with. Jupyter kernels +# launched by an editor do not reliably inherit the variable above, and a +# kernel already running when this script runs keeps its old environment. +# +# The file works out the location itself rather than storing it: .venv lives +# inside the repository, so the parent of sys.prefix is the repository root. +# A stored absolute path would be wrong for everybody else the moment this +# script runs against a different mount point, such as inside a container. +if ! python - <<'PYEOF' +import pathlib +import sys +import sysconfig + +prefix = pathlib.Path(sys.prefix) +if not (prefix / "pyvenv.cfg").exists(): + # Refuse to write anywhere but the tutorial's own environment. Without this + # check a mis-ordered PATH would send the file into the read-only stack. + raise SystemExit(f"{sys.executable} is not the tutorial virtual environment") + +site_packages = pathlib.Path(sysconfig.get_paths()["purelib"]) +site_packages.mkdir(parents=True, exist_ok=True) +(site_packages / "drdcalo-tutorials.pth").write_text( + "import os, site, sys;" + " site.addsitedir(os.path.join(os.path.dirname(sys.prefix), 'python'))\n" +) +PYEOF +then + echo "Warning: could not register drdcalo_tutorials inside .venv." + echo " It stays importable in this shell through PYTHONPATH, but a Jupyter" + echo " kernel started by an editor may not find it. Delete .venv and source" + echo " setup.sh again to rebuild the environment." +fi + +# Step 5 records absolute paths for an editor running on the host. Inside a +# container the repository is mounted somewhere else, so writing .env there +# would point VS Code at paths that do not exist on the host and the notebooks +# would stop working. Everything above applies either way. +_TUTORIAL_IN_CONTAINER=false +if [ -n "${APPTAINER_CONTAINER:-}${SINGULARITY_CONTAINER:-}" ] || \ + [ -n "${APPTAINER_NAME:-}${SINGULARITY_NAME:-}" ] || \ + [ -f /.dockerenv ]; then + _TUTORIAL_IN_CONTAINER=true +fi + +if ${_TUTORIAL_IN_CONTAINER}; then + echo + echo "Container detected: skipping .env." + echo " It is read on the host, and would record ${_TUTORIAL_ROOT}," + echo " which is this container's mount point. Source setup.sh on the host" + echo " to set it up for notebooks opened there." + echo + echo "Environment ready (container)." + echo " python : $(command -v python)" + export DRDCALO_TUTORIALS_ROOT="${_TUTORIAL_ROOT}" + unset _TUTORIAL_ROOT _TUTORIAL_PREV_DIR _TUTORIAL_IN_CONTAINER _TUTORIAL_NEED_STACK + return 0 +fi + +# 5. Record the part of the environment the notebooks need, so that Jupyter +# kernels started by an editor, which do not inherit this shell, still find +# ROOT, DD4hep and podio. The repository location is written down too, so that +# notebooks can locate their input file whatever directory the editor happens +# to start them in. +export DRDCALO_TUTORIALS_ROOT="${_TUTORIAL_ROOT}" +if ! DRDCALO_ENV_FILE="${_TUTORIAL_ROOT}/.env" python - <<'PYEOF' +import os +import stat + +# Only the variables that make the key4hep runtime work are written out. The +# shell environment as a whole is not ours to copy into a file: it routinely +# carries Kerberos ticket paths, proxy credentials, CI tokens and similar, and +# .env is read by an editor rather than by a login shell. Everything else the +# kernel needs (HOME, USER, ...) it already inherits from the editor. +# +# If a notebook ever reports something missing from its environment, add the +# variable here rather than widening this into "copy everything". +KEEP = { + # Executables, shared libraries and Python packages + "PATH", + "LD_LIBRARY_PATH", + "PYTHONPATH", + "CMAKE_PREFIX_PATH", + # This repository: the virtual environment and the input-file resolution + "VIRTUAL_ENV", + "VIRTUAL_ENV_PROMPT", + "DRDCALO_TUTORIALS_ROOT", + # ROOT and its Python bindings, which is what `import dd4hep` goes through + "ROOTSYS", + "ROOT_INCLUDE_PATH", + "ROOT_LIBRARY_PATH", + "ROOT_VERSION", + "CLING_STANDARD_PCH", + "CPPYY_API_PATH", + "CPPYY_BACKEND_LIBRARY", + # The stack itself + "KEY4HEP_STACK", + "key4hep_stack_version", + # Plugins, geometry and detector data used by the exercises + "GAUDI_PLUGIN_PATH", + "MARLIN_DLL", + "DD4HEP", + "DD4hepINSTALL", + "DD4hep_DIR", + "DD4hep_ROOT", + "K4GEO", + "k4geo_DIR", + "LCGEO", + "lcgeo_DIR", + "LCIO", + # Jupyter, so that a kernel started from this environment finds its own data + "JUPYTER_PATH", + "JUPYTER_CONFIG_PATH", + "JUPYTERLAB_DIR", +} + +# Geant4 dataset locations: G4LEDATA, G4ENSDFSTATEDATA and a dozen more, all +# named the same way and all pointing into the stack. +KEEP_PREFIXES = ("G4",) + +lines = [] +for name, value in sorted(os.environ.items()): + if name not in KEEP and not name.startswith(KEEP_PREFIXES): + continue + if "\n" in value: # cannot be represented in a .env file + continue + lines.append(f'{name}="{value}"\n') + +# Written by hand rather than through a shell redirection so that the file is +# created 0600 instead of whatever the umask happens to be: it describes this +# user's stack and belongs to this user. +path = os.environ["DRDCALO_ENV_FILE"] +fd = os.open(path, os.O_WRONLY | os.O_CREAT | os.O_TRUNC, 0o600) +with os.fdopen(fd, "w") as env_file: + env_file.writelines(lines) +# O_CREAT leaves the mode of an already existing file alone, so set it explicitly. +os.chmod(path, stat.S_IRUSR | stat.S_IWUSR) +PYEOF +then + echo "Failed to write ${_TUTORIAL_ROOT}/.env." + return 1 +fi + +echo +echo "Environment ready." +echo " python : $(command -v python)" +echo " .env : ${_TUTORIAL_ROOT}/.env ($(wc -l < "${_TUTORIAL_ROOT}/.env") variables)" +echo +echo "In VS Code, open ${_TUTORIAL_ROOT} as the workspace folder and select" +echo "the .venv interpreter. VS Code applies .env from the workspace folder," +echo "which is what puts the key4hep runtime into the notebook." + +unset _TUTORIAL_ROOT _TUTORIAL_PREV_DIR _TUTORIAL_IN_CONTAINER _TUTORIAL_NEED_STACK diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 2c3200b..daf6335 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -17,17 +17,35 @@ See the License for the specific language governing permissions and limitations under the License. ]] +set(SIM_FILE ${PROJECT_SOURCE_DIR}/DD4hepTutorials/data/simplecalo2_sample.root) +if(NOT EXISTS ${SIM_FILE}) + message(FATAL_ERROR "Bundled tutorial input not found: ${SIM_FILE}") +endif() + +# Everything the tests write goes into the build directory: the option files are +# named absolutely and the working directory is the build tree, so a `ctest` run +# leaves the checkout untouched. +set(GAUDI_OPTIONS ${PROJECT_SOURCE_DIR}/GaudiTutorial) + add_test(NAME test_EventStats - WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} - COMMAND k4run GaudiTutorial/EventStats/options/runEventStatsSolution.py --IOSvc.Input ${PROJECT_SOURCE_DIR}/GaudiTutorial/data/simpleCalo_simulation.root --IOSvc.Output ${PROJECT_SOURCE_DIR}/GaudiTutorial/data/simpleCalo_eventStats.root) + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMAND k4run ${GAUDI_OPTIONS}/EventStats/options/runEventStatsSolution.py --IOSvc.Input ${SIM_FILE} --IOSvc.Output ${CMAKE_CURRENT_BINARY_DIR}/simpleCalo_eventStats.root --EventStats.HistogramFile ${CMAKE_CURRENT_BINARY_DIR}/debugEventStats.root) + +# The multithreaded configuration is what motivates the mutex around the shared +# histograms in EventStats.cpp, so it is exercised as well. Ten events across six +# threads is cheap, and ThreadSanitizer aside it is what catches an unprotected +# TH1::Fill turning into a crash or a corrupted histogram. +add_test(NAME test_EventStatsMultithreading + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMAND k4run ${GAUDI_OPTIONS}/EventStats/options/runEventStatsMultithreading.py --IOSvc.Input ${SIM_FILE} --IOSvc.Output ${CMAKE_CURRENT_BINARY_DIR}/simpleCalo_eventStats_multithreaded.root --EventStats.HistogramFile ${CMAKE_CURRENT_BINARY_DIR}/debugEventStats_multithreaded.root) add_test(NAME test_MoliereRadius - WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} - COMMAND k4run GaudiTutorial/MoliereRadiusFunctional/options/runMoliereRadiusSolution.py --IOSvc.Input ${PROJECT_SOURCE_DIR}/GaudiTutorial/data/simpleCalo_simulation.root --IOSvc.Output ${PROJECT_SOURCE_DIR}/GaudiTutorial/data/simpleCalo_moliereRadius.root) + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMAND k4run ${GAUDI_OPTIONS}/MoliereRadiusFunctional/options/runMoliereRadiusSolution.py --IOSvc.Input ${SIM_FILE} --IOSvc.Output ${CMAKE_CURRENT_BINARY_DIR}/simpleCalo_moliereRadius.root) add_test(NAME test_RandomNoiseDigitizer - WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} - COMMAND k4run GaudiTutorial/RandomNoiseDigitizer/options/runRandomNoiseDigitizerSolution.py --IOSvc.Input ${PROJECT_SOURCE_DIR}/GaudiTutorial/data/simpleCalo_simulation.root --IOSvc.Output ${PROJECT_SOURCE_DIR}/GaudiTutorial/data/simpleCalo_noiseDigitizer.root) + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMAND k4run ${GAUDI_OPTIONS}/RandomNoiseDigitizer/options/runRandomNoiseDigitizerSolution.py --IOSvc.Input ${SIM_FILE} --IOSvc.Output ${CMAKE_CURRENT_BINARY_DIR}/simpleCalo_noiseDigitizer.root) # For LD_LIBRARY_PATH # ${PROJECT_BINARY_DIR} # The files .components and .confdb are here @@ -46,20 +64,35 @@ add_test(NAME test_RandomNoiseDigitizer # https://github.com/key4hep/k4FWCore/blob/42497b71c23425a3a26648f30a1285a499dcda28/test/k4FWCoreTest/CMakeLists.txt#L71 -execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory - ${PROJECT_SOURCE_DIR}/GaudiTutorial/data) - -set(SIM_FILE ${PROJECT_SOURCE_DIR}/GaudiTutorial/data/simpleCalo_simulation.root) +# Check that every notebook is well formed and that its code cells are valid Python. +add_test(NAME test_NotebooksValid + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMAND python ${CMAKE_CURRENT_SOURCE_DIR}/check_notebooks.py) -if(NOT EXISTS ${SIM_FILE}) - message(STATUS "Downloading simpleCalo_simulation.root ...") - execute_process( - COMMAND wget -O ${SIM_FILE} https://fccsw.web.cern.ch/fccsw/tutorials/DRD6_2025_Ancona/simpleCalo_simulation.root - ) -else() - message(STATUS "File already exists: ${SIM_FILE}, skipping download.") +# Execute the solution notebook against the bundled sample by default. Developers can still +# override SIMPLECALO2_FILE at configure time to exercise a full simulation output. +if(NOT DEFINED SIMPLECALO2_FILE) + set(SIMPLECALO2_FILE + ${PROJECT_SOURCE_DIR}/DD4hepTutorials/data/simplecalo2_sample.root) +endif() +if(NOT EXISTS ${SIMPLECALO2_FILE}) + message(FATAL_ERROR "EDM4hep notebook input not found: ${SIMPLECALO2_FILE}") endif() +add_test(NAME test_ReadEdm4hepNotebook + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMAND jupyter nbconvert --to notebook --execute + ${PROJECT_SOURCE_DIR}/DD4hepTutorials/notebooks/readEdm4hepSolution.ipynb + --output ${CMAKE_CURRENT_BINARY_DIR}/readEdm4hepSolution_executed.ipynb) + +# The Section 1 notebook runs on the same bundled file, which carries the same +# simplecaloRO collection. SIMPLECALO1_FILE keeps it off any full simulation that +# happens to be sitting in the source tree. +add_test(NAME test_CellEnergySumNotebook + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMAND jupyter nbconvert --to notebook --execute + ${PROJECT_SOURCE_DIR}/DD4hepTutorials/notebooks/plot_cell_energy_sum.ipynb + --output ${CMAKE_CURRENT_BINARY_DIR}/plot_cell_energy_sum_executed.ipynb) get_filename_component(CURRENT_DIR_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME) @@ -75,6 +108,7 @@ ${CMAKE_CURRENT_BINARY_DIR}/genConfDir/${CURRENT_DIR_NAME}:\ $ENV{LD_LIBRARY_PATH};\ \ PYTHONPATH=\ +${PROJECT_SOURCE_DIR}/python:\ ${PROJECT_BINARY_DIR}/GaudiTutorial/EventStats/genConfDir:\ ${PROJECT_BINARY_DIR}/GaudiTutorial/MoliereRadiusFunctional/genConfDir:\ ${PROJECT_BINARY_DIR}/GaudiTutorial/RandomNoiseDigitizer/genConfDir:\ @@ -85,6 +119,16 @@ $ENV{PYTHONPATH};\ get_property(test_names DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} PROPERTY TESTS) set_tests_properties(${test_names} PROPERTIES ENVIRONMENT "${test_environment}") +# Must come after set_tests_properties() above, which would otherwise overwrite it. +if(TEST test_ReadEdm4hepNotebook) + set_property(TEST test_ReadEdm4hepNotebook APPEND PROPERTY + ENVIRONMENT "SIMPLECALO2_FILE=${SIMPLECALO2_FILE}") +endif() +if(TEST test_CellEnergySumNotebook) + set_property(TEST test_CellEnergySumNotebook APPEND PROPERTY + ENVIRONMENT "SIMPLECALO1_FILE=${SIM_FILE}") +endif() + message(STATUS "================ Test Environment ================") message(STATUS "${test_environment}") -message(STATUS "=================================================") \ No newline at end of file +message(STATUS "=================================================") diff --git a/test/check_notebooks.py b/test/check_notebooks.py new file mode 100755 index 0000000..b51a619 --- /dev/null +++ b/test/check_notebooks.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python3 +# +# Copyright (c) 2020-2024 Key4hep-Project. +# +# This file is part of Key4hep. +# See https://key4hep.github.io/key4hep-doc/ for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +"""Validate the tutorial notebooks. + +Checks that every notebook is well formed and that all of its code cells are +syntactically valid Python. This runs everywhere, including where no simulation +output is available to execute the notebooks against. +""" +import sys +from pathlib import Path + +import nbformat + + +def check(path): + problems = [] + nb = nbformat.read(path, as_version=4) + nbformat.validate(nb) + + n_code = 0 + for number, cell in enumerate(nb.cells): + if cell.cell_type != "code": + continue + n_code += 1 + # IPython magics (%jsroot on) and shell escapes are not valid Python + source = "\n".join( + "" if line.lstrip().startswith(("%", "!")) else line + for line in cell.source.splitlines() + ) + try: + compile(source, f"{path.name}:cell{number}", "exec") + except SyntaxError as exc: + problems.append(f"{path.name} cell {number}: {exc}") + + if n_code == 0: + problems.append(f"{path.name}: contains no code cells") + return n_code, problems + + +def main(): + root = Path(__file__).resolve().parent.parent + notebooks = sorted(root.glob("*/notebooks/*.ipynb")) + if not notebooks: + print("ERROR: no notebooks found", file=sys.stderr) + return 1 + + failed = [] + for path in notebooks: + n_code, problems = check(path) + status = "OK" if not problems else "FAIL" + print(f"[{status}] {path.relative_to(root)} ({n_code} code cells)") + failed.extend(problems) + + for problem in failed: + print(" " + problem, file=sys.stderr) + return 1 if failed else 0 + + +if __name__ == "__main__": + sys.exit(main())