Add Bethe Bloch energy loss model (ATTPC#259)

* Add CLAUDE.md file

* Add AtELossBetheBloch analytic energy loss model

Implements the Bethe-Bloch stopping power formula as a new AtELossModel
subclass supporting heavy charged particles (PDG 2022, Eq. 34.1) and
electrons (Leo 1994, Eq. 2.38). A tk::spline cache of dx/dE enables
O(log n) range integrals. Newton's method is used for GetEnergy,
matching the AtELossTable strategy. Straggling uses the Bohr
approximation. Low-energy breakdown of the formula is handled by
extrapolating dEdx ∝ sqrt(E) below the validity threshold.

Includes 8 unit tests (proton/alpha/pion/electron in H2 and Ar,
self-consistency, straggling sanity, Bloch approximation) and a ROOT
macro to plot proton stopping power in Al for visual comparison against
the PDG reference curve.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

* Add more robust pion tests

* Fix AtELossBetheBloch setters to rebuild spline automatically

SetMaterial, SetI, and SetDensity all change physical parameters that
the cached dx/dE spline depends on, but previously left it stale.
SetMaterial now calls BuildSpline at the end (removing the redundant
explicit call from the constructor). SetI is moved from an inline
header one-liner to a proper definition that rebuilds the spline.
SetDensity overrides the base class virtual and rebuilds the spline
after updating fDensity.

Adds SettersRebuildSpline test to verify each setter takes effect
immediately without a manual BuildSpline call.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

* Add record of Claude log for Bethe Impl

* Move level of debug statement

* Fix clang-formatting

---------

Co-authored-by: Claude Sonnet 4.6 <noreply@anthropic.com>

Author: anthoak13

.claude/2026-03-03-131255-implement-beth-bloch.txt

@@ -0,0 +1,146 @@
+# CLAUDE.md
+
+This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
+
+## About
+
+ATTPCROOT is a ROOT/FairRoot-based C++ framework for simulation and analysis of Active Target Time Projection Chamber (AT-TPC) detector data. It integrates with FairSoft (provides ROOT, Geant4, VMC) and FairRoot (provides the task/run framework).
+
+## Environment Setup
+
+Before building or running anything, the environment must be loaded. The VSCode terminal auto-sources this on startup:
+
+```bash
+source build/config.sh
+```
+
+This sets `LD_LIBRARY_PATH`, `ROOTSYS`, `VMCWORKDIR`, `ROOT_INCLUDE_PATH`, and Geant4 data paths. The key install paths on this machine are:
+- FairRoot: `~/fair_install/FairRootInstall`
+- FairSoft: `~/fair_install/FairSoftInstall`
+
+For CMake configuration (not the build itself), the following env vars are needed:
+```bash
+export FAIRROOTPATH=~/fair_install/FairRootInstall
+export SIMPATH=~/fair_install/FairSoftInstall
+```
+
+## Build Commands
+
+```bash
+# Configure (from repo root, out-of-source into build/)
+cmake -S . -B build -DCMAKE_PREFIX_PATH=~/fair_install/hdf5
+
+# Build (10 parallel jobs)
+cmake --build build -j10
+
+# Build a specific target
+cmake --build build --target AtSimulationData -j10
+
+# Run all unit tests
+cd build && ctest -V
+
+# Run a specific test binary directly
+./build/tests/AtSimulationDataTests
+./build/tests/AtGeneratorsTests
+./build/tests/AtToolsTests
+```
+
+Tests are built by default (`BUILD_TESTS=ON`). Test binaries are placed in `build/tests/`.
+
+## Code Formatting
+
+The project uses `clang-format-17` with the config in `.clang-format` (based on LLVM style, 3-space indent, 120-char column limit). Format on save is configured in VSCode. To format manually:
+
+```bash
+clang-format-17 -i <file>
+# Or format all changed files:
+scripts/formatAll.sh
+```
+
+Static analysis uses `clang-tidy` (config in `.clang-tidy`): modernize-* and cppcoreguidelines-* checks.
+
+## Architecture
+
+The framework follows FairRoot's task pipeline pattern: data flows through a chain of `FairTask` subclasses, persisted as ROOT TClonesArrays in a TTree.
+
+### Module Overview
+
+| Module | Purpose |
+|--------|---------|
+| `AtSimulationData` | MC truth data: `AtStack`, `AtMCTrack`, `AtMCPoint`, `AtVertexPropagator` |
+| `AtGenerators` | FairGenerator subclasses for beam/reaction/decay event generation |
+| `AtDetectors` | Geant4 sensitive detector implementations (AT-TPC, GADGET-II, SpecMAT, etc.) |
+| `AtDigitization` | Converts MC points → simulated pad signals (`AtClusterize`, `AtPulse`) |
+| `AtUnpack` | Unpacks raw experimental data (GET/GRAW, HDF5, ROOT formats) |
+| `AtData` | Core data classes: `AtRawEvent`, `AtEvent`, `AtHit`, `AtPad`, `AtTrack` |
+| `AtMap` | Pad mapping between electronics channels and detector geometry |
+| `AtParameter` | FairRuntimeDb parameter containers |
+| `AtReconstruction` | PSA, pattern recognition, fitting tasks |
+| `AtTools` | Utilities: energy loss (CATIMA), kinematics, space charge, hit sampling |
+| `AtAnalysis` | High-level analysis tasks and `AtRunAna` |
+| `AtS800` | S800 spectrograph data handling |
+| `AtEventDisplay` | ROOT Eve-based event display |
+
+### Data Flow
+
+**Simulation pipeline** (macros in `macro/Simulation/`):
+1. `FairPrimaryGenerator` with an `AtReactionGenerator` subclass → particle stack
+2. Geant4/VMC transport → `AtMCPoint` hits in sensitive volumes
+3. `AtClusterizeTask` → electron clusters from ionization
+4. `AtPulseTask` → simulated pad signals (`AtRawEvent`)
+
+**Reconstruction pipeline** (macros in `macro/Unpack_*/` and `macro/Analysis/`):
+1. `AtUnpackTask` → `AtRawEvent` (raw pad traces)
+2. `AtFilterTask` → filtered `AtRawEvent`
+3. `AtPSAtask` (Pulse Shape Analysis) → `AtEvent` with `AtHit` objects
+4. `AtSampleConsensusTask` / `AtPRAtask` → `AtPatternEvent` with `AtTrack` objects
+5. `AtMCFitterTask` / `AtFitterTask` → fitted tracks with kinematics
+
+### Key Design Patterns
+
+**FairTask subclasses**: Each processing step is a `FairTask`. They retrieve input branches via `TClonesArray*` from `FairRootManager` in `Init()` and process them in `Exec()`.
+
+**AtReactionGenerator**: Abstract base for all reaction generators. Subclasses implement `GenerateReaction()`. The `ReadEvent()` method is `final` and handles the beam/reaction event alternation via `AtVertexPropagator`. Generators can be chained for sequential decays using `SetSequentialDecay(true)`.
+
+**AtVertexPropagator**: Singleton (`AtVertexPropagator::Instance()`) that communicates vertex, momentum, and kinematics between chained generators and downstream tasks. Alternates beam/reaction events via `EndEvent()`. Use `ResetForTesting()` in unit tests.
+
+**PSA (Pulse Shape Analysis)**: `AtPSA` is the abstract base. Concrete implementations (`AtPSAMax`, `AtPSAFull`, `AtPSADeconv`, etc.) extract hits from pad traces. `AtPSAComposite` allows chaining PSA methods.
+
+**Energy loss**: `AtELossModel` is the base; `AtELossCATIMA` wraps the CATIMA library (fetched automatically if not installed). Used via `AtELossManager`.
+
+### Adding a New Module Library
+
+Each module's `CMakeLists.txt` follows this pattern:
+```cmake
+set(LIBRARY_NAME MyModule)
+set(SRCS file1.cxx file2.cxx)
+set(DEPENDENCIES ATTPCROOT::AtData ROOT::Core ...)
+set(TEST_SRCS MyModuleTest.cxx)
+attpcroot_generate_tests(${LIBRARY_NAME}Tests SRCS ${TEST_SRCS} DEPS ${LIBRARY_NAME})
+generate_target_and_root_library(${LIBRARY_NAME} LINKDEF ${LIBRARY_NAME}LinkDef.h SRCS ${SRCS} DEPS_PUBLIC ${DEPENDENCIES})
+```
+
+ROOT dictionary generation requires a `*LinkDef.h` file listing every class that needs ROOT reflection. Every class with `ClassDef` in its header must appear in the LinkDef.
+
+**LinkDef streamer suffix** — the suffix after the class name controls whether ROOT generates a full I/O streamer:
+
+- `ClassName +;` — generates a full streamer. Use this **only** for classes that are written to a ROOT file (i.e. stored as a branch in a `TClonesArray` or `TTree`). Examples: `AtHit`, `AtEvent`, `AtRawEvent`, `AtMCPoint`.
+- `ClassName -!;` — registers the class for reflection (usable in interpreted macros, usable as a pointer type) but **suppresses the streamer**. Use this for every class that is never written to disk: tasks, algorithms, models, samplers, etc. Examples: `AtELossModel`, `AtSpaceChargeModel`, `AtSample` subclasses.
+
+Generating an unnecessary streamer bloats the dictionary and binary with unused I/O code, so the default should be `-!` unless disk persistence is actually required.
+
+**C++ vs ROOT typedefs** — use plain C++ types (`bool`, `int`, `double`, `std::string`) in all classes that are not written to a ROOT file. Only use ROOT's portability typedefs (`Bool_t`, `Int_t`, `Double_t`, etc.) in classes that are persisted to disk, where ROOT's cross-platform I/O guarantees matter. Mixing ROOT types into purely algorithmic or task classes is unnecessary overhead.
+
+### Unit Tests
+
+Tests use Google Test, placed alongside source files (e.g., `AtVertexPropagatorTest.cxx` in `AtSimulationData/`). Register them in the module's `CMakeLists.txt` via `attpcroot_generate_tests()`. Tests must not access external files or network resources.
+
+### Macros
+
+ROOT macros (`.C` files) in `macro/` are the primary user interface for running simulations and analyses. They are not compiled into libraries—they run interpreted by ROOT. Integration tests live in `macro/tests/`.
+
+## Contributing
+
+- PRs must target the `develop` branch and apply cleanly (fast-forward, no merge commits).
+- Commit messages: present imperative mood, ≤72 characters summary.
+- All PRs are checked by `clang-format`, `clang-tidy`, and the unit test suite.

.claude/CLAUDE.md

@@ -0,0 +1,221 @@
+#include "AtELossBetheBloch.h"
+
+#include <FairLogger.h>
+
+#include <cmath>
+#include <vector>
+
+namespace AtTools {
+
+AtELossBetheBloch::AtELossBetheBloch(double part_q, double part_mass, int mat_Z, int mat_A, double density, double I_eV)
+   : AtELossModel(density), fPart_q(part_q), fPart_mass(part_mass)
+{
+   SetMaterial(mat_Z, mat_A, density, I_eV);
+}
+
+void AtELossBetheBloch::SetMaterial(int mat_Z, int mat_A, double density, double I_eV)
+{
+   fMat_Z = mat_Z;
+   fMat_A = mat_A;
+   fDensity = density;
+   if (I_eV <= 0)
+      fI_MeV = 13.5 * mat_Z * 1e-6; // Bloch approximation: I ≈ 13.5·Z eV
+   else
+      fI_MeV = I_eV * 1e-6;
+   BuildSpline();
+}
+
+void AtELossBetheBloch::SetI(double I_eV)
+{
+   fI_MeV = I_eV * 1e-6;
+   BuildSpline();
+}
+
+void AtELossBetheBloch::SetDensity(double density)
+{
+   fDensity = density;
+   BuildSpline();
+}
+
+void AtELossBetheBloch::BuildSpline(double E_min_MeV, double E_max_MeV, int nPoints)
+{
+   std::vector<double> energies(nPoints);
+   std::vector<double> dedxValues(nPoints);
+
+   double logMin = std::log(E_min_MeV);
+   double logMax = std::log(E_max_MeV);
+
+   for (int i = 0; i < nPoints; ++i) {
+      double t = static_cast<double>(i) / (nPoints - 1);
+      energies[i] = std::exp(logMin + t * (logMax - logMin));
+      dedxValues[i] = GetdEdx_formula(energies[i]);
+   }
+
+   // The Bethe-Bloch formula breaks down at low energies where the log argument drops below 1.
+   // Find the lowest-energy valid grid point, then extrapolate below it using dEdx ∝ sqrt(E),
+   // consistent with the low-energy Lindhard stopping power behavior.
+   int firstValid = -1;
+   for (int i = 0; i < nPoints; ++i) {
+      if (dedxValues[i] > 0) {
+         firstValid = i;
+         break;
+      }
+   }
+
+   if (firstValid < 0) {
+      LOG(error) << "Bethe-Bloch formula returned zero for all energy grid points!";
+      return;
+   }
+   if (firstValid > 0) {
+      LOG(debug) << "Bethe-Bloch formula invalid below " << energies[firstValid]
+                 << " MeV; extrapolating with dEdx ∝ sqrt(E).";
+      for (int i = 0; i < firstValid; ++i)
+         dedxValues[i] = dedxValues[firstValid] * std::sqrt(energies[i] / energies[firstValid]);
+   }
+
+   std::vector<double> dXdE(nPoints);
+   for (int i = 0; i < nPoints; ++i)
+      dXdE[i] = (dedxValues[i] > 0) ? 1.0 / dedxValues[i] : 1e10;
+
+   fdXdE = tk::spline(energies, dXdE);
+}
+
+double AtELossBetheBloch::GetdEdx_formula(double energy) const
+{
+   if (IsElectron())
+      return GetdEdx_electron(energy);
+   return GetdEdx_heavy(energy);
+}
+
+double AtELossBetheBloch::GetdEdx_heavy(double energy) const
+{
+   if (energy <= 0)
+      return 0;
+
+   double M = fPart_mass;
+   double gamma = 1.0 + energy / M;
+   double beta2 = 1.0 - 1.0 / (gamma * gamma);
+
+   if (beta2 <= 0)
+      return 0;
+
+   // Maximum kinetic energy transferable in a single collision (PDG 2022, Eq. 34.2)
+   double Tmax = 2.0 * kM_e * beta2 * gamma * gamma / (1.0 + 2.0 * gamma * kM_e / M + (kM_e / M) * (kM_e / M));
+
+   // Argument of logarithm: 2mₑc²β²γ²Tmax / I²
+   double logArg = 2.0 * kM_e * beta2 * gamma * gamma * Tmax / (fI_MeV * fI_MeV);
+   if (logArg <= 0)
+      return 0;
+
+   double z = fPart_q;
+   // Standard Bethe-Bloch in MeV/cm (PDG 2022, Eq. 34.1), density correction neglected
+   double dedx_cm =
+      kK * z * z * (static_cast<double>(fMat_Z) / fMat_A) * fDensity / beta2 * (0.5 * std::log(logArg) - beta2);
+
+   if (dedx_cm <= 0)
+      return 0;
+
+   return dedx_cm / 10.0; // Convert MeV/cm → MeV/mm
+}
+
+double AtELossBetheBloch::GetdEdx_electron(double energy) const
+{
+   if (energy <= 0)
+      return 0;
+
+   double tau = energy / kM_e;
+   double beta2 = tau * (tau + 2.0) / ((tau + 1.0) * (tau + 1.0));
+
+   if (beta2 <= 0)
+      return 0;
+
+   // Møller exchange correction (Leo 1994, Eq. 2.38)
+   double Fminus = (1.0 + tau * tau / 8.0 - (2.0 * tau + 1.0) * std::log(2.0)) / ((tau + 1.0) * (tau + 1.0));
+
+   // Argument of logarithm: τ√(τ+2)·mₑc² / (√2·I)
+   double logArg = tau * std::sqrt(tau + 2.0) * kM_e / (std::sqrt(2.0) * fI_MeV);
+   if (logArg <= 0)
+      return 0;
+
+   // Modified Bethe formula for electrons in MeV/cm
+   double dedx_cm = kK * (static_cast<double>(fMat_Z) / fMat_A) * fDensity / beta2 * (std::log(logArg) + Fminus);
+
+   if (dedx_cm <= 0)
+      return 0;
+
+   return dedx_cm / 10.0; // Convert MeV/cm → MeV/mm
+}
+
+double AtELossBetheBloch::GetdEdx(double energy) const
+{
+   return 1.0 / fdXdE(energy);
+}
+
+double AtELossBetheBloch::GetRange(double energyIni, double energyFin) const
+{
+   if (energyIni == energyFin)
+      return 0;
+   if (energyFin < fdXdE.get_x_min()) {
+      LOG(debug) << "Attempting to integrate energy to " << energyFin << " when min energy in table is "
+                 << fdXdE.get_x_min();
+      energyFin = fdXdE.get_x_min();
+   }
+   return fdXdE.integrate(energyFin, energyIni);
+}
+
+double AtELossBetheBloch::GetEnergyLoss(double energyIni, double distance) const
+{
+   return energyIni - GetEnergy(energyIni, distance);
+}
+
+double AtELossBetheBloch::GetEnergy(double energyIni, double distance) const
+{
+   if (distance == 0)
+      return energyIni;
+   if (energyIni < 1e-6 || GetRange(energyIni) < distance)
+      return 0.0;
+
+   const int maxIt = 100;
+   const double distErr = 1e-4; // mm
+
+   double guessEnergy = energyIni - GetdEdx(energyIni) * distance;
+   for (int i = 0; i < maxIt; ++i) {
+      double range = GetRange(energyIni, guessEnergy);
+      if (std::fabs(range - distance) < distErr) {
+         LOG(debug) << "Energy converged in " << i + 1 << " iterations.";
+         return guessEnergy;
+      }
+      guessEnergy += GetdEdx(guessEnergy) * (range - distance);
+      if (guessEnergy < fdXdE.get_x_min() * 1.01) {
+         guessEnergy = fdXdE.get_x_min();
+         return guessEnergy;
+      }
+   }
+
+   LOG(error) << "Energy calculation (" << energyIni << " MeV through " << distance << " mm) failed to converge in "
+              << maxIt << " iterations!";
+   return -1;
+}
+
+double AtELossBetheBloch::GetElossStraggling(double energyIni, double energyFin) const
+{
+   double dx_mm = GetRange(energyIni, energyFin);
+   if (dx_mm <= 0)
+      return 0;
+
+   // Bohr approximation: σ² = K·z²·(Z/A)·ρ·Δx[cm]·mₑc²
+   double z = IsElectron() ? 1.0 : fPart_q;
+   double dx_cm = dx_mm * 0.1;
+   double sigma2 = kK * z * z * (static_cast<double>(fMat_Z) / fMat_A) * fDensity * dx_cm * kM_e;
+   return std::sqrt(sigma2);
+}
+
+double AtELossBetheBloch::GetdEdxStraggling(double energyIni, double energyFin) const
+{
+   double dx_mm = GetRange(energyIni, energyFin);
+   if (dx_mm <= 0)
+      return 0;
+   return GetElossStraggling(energyIni, energyFin) / dx_mm;
+}
+
+} // namespace AtTools