diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index 6e3f1de01..138117faa 100644 --- a/.claude/CLAUDE.md +++ b/.claude/CLAUDE.md @@ -1,146 +1,60 @@ # 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: +ATTPCROOT is a ROOT/FairRoot-based C++ framework for simulation and analysis of Active Target Time Projection Chamber (AT-TPC) detector data. + +## Documentation + +Full developer documentation lives in `docs/`. See [docs/index.md](../docs/index.md) for the full map. Quick topic links: + +| Topic | File | +|-------|------| +| First-time install | [tooling/installation.md](../docs/tooling/installation.md) | +| Daily use (build/test) | [tooling/daily-use.md](../docs/tooling/daily-use.md) | +| Testing patterns | [tooling/testing.md](../docs/tooling/testing.md) | +| Contributor guide | [contributing/guide.md](../docs/contributing/guide.md) | +| Adding a new module | [contributing/new-module.md](../docs/contributing/new-module.md) | +| Code style | [contributing/code-style.md](../docs/contributing/code-style.md) | +| Module overview | [reference/modules.md](../docs/reference/modules.md) | +| Data model | [reference/data-model.md](../docs/reference/data-model.md) | +| Branch I/O contracts | [reference/branch-io-contracts.md](../docs/reference/branch-io-contracts.md) | +| Simulation pipeline | [subsystems/simulation-pipeline.md](../docs/subsystems/simulation-pipeline.md) | +| Reconstruction pipeline | [subsystems/reconstruction-pipeline.md](../docs/subsystems/reconstruction-pipeline.md) | +| Event generators | [subsystems/generators.md](../docs/subsystems/generators.md) | +| Pulse shape analysis | [subsystems/psa.md](../docs/subsystems/psa.md) | +| Energy loss | [subsystems/energy-loss.md](../docs/subsystems/energy-loss.md) | + +## Quick Reference: Build & Test ```bash -clang-format-17 -i -# Or format all changed files: -scripts/formatAll.sh +source build/config.sh # load environment (do this first) +cmake --build build -j10 # build everything +cd build && ctest -V # run all unit tests ``` -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}) -``` +## Code-Writing Rules -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. +These rules apply whenever editing or adding C++ code in this repo. -**LinkDef streamer suffix** — the suffix after the class name controls whether ROOT generates a full I/O streamer: +### LinkDef Streamer Suffixes -- `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. +Every class in a `*LinkDef.h` file must use the correct suffix: -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. +- `ClassName +;` — generates a full I/O streamer. **Only** for classes written to a ROOT file (stored in a `TClonesArray` or `TTree` branch). Examples: `AtHit`, `AtEvent`, `AtRawEvent`, `AtMCPoint`. +- `ClassName -!;` — reflection only, no streamer. Use for **everything else**: tasks, algorithms, models, samplers. Examples: `AtPSAMax`, `AtELossModel`, `AtFitterTask`. -**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. +Default to `-!` unless disk persistence is actually required. -### Unit Tests +### C++ vs ROOT Typedefs -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. +- Non-persisted classes (tasks, algorithms, models): use `bool`, `int`, `double`, `std::string`. +- Persisted data classes only: use `Bool_t`, `Int_t`, `Double_t`, etc. -### Macros +### Test Isolation -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/`. +Unit tests must not access external files or network resources. Hardcode test data inline. ## 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. +- PRs target the `develop` branch; fast-forward only (no merge commits). +- Commit messages: present imperative mood, ≤72 characters. +- All PRs must pass `clang-format`, `clang-tidy`, and unit tests. diff --git a/.github/workflows/CI-tests.yml b/.github/workflows/CI-tests.yml index 810be15de..33b8a9e01 100644 --- a/.github/workflows/CI-tests.yml +++ b/.github/workflows/CI-tests.yml @@ -18,7 +18,7 @@ jobs: runs-on: ubuntu-latest container: anthoak13/attpcroot-deps steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - name: Configure CMake # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index ffb510443..56253261b 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -8,7 +8,7 @@ jobs: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v6 - name: Run clang-format uses: DoozyX/clang-format-lint-action@v0.17 # These are optional (defaults displayed) diff --git a/.gitignore b/.gitignore index f0c9c9839..081c538ef 100755 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ +*.pdf *.swp *.root *.log @@ -57,3 +58,6 @@ macros/data .DS_Store .nfs* compiled/E20009Analysis/.NabinWorkDir/ + +.claude/plans/* + diff --git a/AGENTS.md b/AGENTS.md new file mode 100644 index 000000000..3b59135d5 --- /dev/null +++ b/AGENTS.md @@ -0,0 +1,60 @@ +# AGENTS.md + +ATTPCROOT is a ROOT/FairRoot-based C++ framework for simulation and analysis of Active Target Time Projection Chamber (AT-TPC) detector data. + +## Documentation + +Full developer documentation lives in `docs/`. See [docs/index.md](docs/index.md) for the full map. Quick topic links: + +| Topic | File | +|-------|------| +| First-time install | [tooling/installation.md](docs/tooling/installation.md) | +| Daily use (build/test) | [tooling/daily-use.md](docs/tooling/daily-use.md) | +| Testing patterns | [tooling/testing.md](docs/tooling/testing.md) | +| Contributor guide | [contributing/guide.md](docs/contributing/guide.md) | +| Adding a new module | [contributing/new-module.md](docs/contributing/new-module.md) | +| Code style | [contributing/code-style.md](docs/contributing/code-style.md) | +| Module overview | [reference/modules.md](docs/reference/modules.md) | +| Data model | [reference/data-model.md](docs/reference/data-model.md) | +| Branch I/O contracts | [reference/branch-io-contracts.md](docs/reference/branch-io-contracts.md) | +| Simulation pipeline | [subsystems/simulation-pipeline.md](docs/subsystems/simulation-pipeline.md) | +| Reconstruction pipeline | [subsystems/reconstruction-pipeline.md](docs/subsystems/reconstruction-pipeline.md) | +| Event generators | [subsystems/generators.md](docs/subsystems/generators.md) | +| Pulse shape analysis | [subsystems/psa.md](docs/subsystems/psa.md) | +| Energy loss | [subsystems/energy-loss.md](docs/subsystems/energy-loss.md) | + +## Quick Reference: Build & Test + +```bash +source build/config.sh # load environment (do this first) +cmake --build build -j10 # build everything +cd build && ctest -V # run all unit tests +``` + +## Code-Writing Rules + +These rules apply whenever editing or adding C++ code in this repo. + +### LinkDef Streamer Suffixes + +Every class in a `*LinkDef.h` file must use the correct suffix: + +- `ClassName +;` — generates a full I/O streamer. **Only** for classes written to a ROOT file (stored in a `TClonesArray` or `TTree` branch). Examples: `AtHit`, `AtEvent`, `AtRawEvent`, `AtMCPoint`. +- `ClassName -!;` — reflection only, no streamer. Use for **everything else**: tasks, algorithms, models, samplers. Examples: `AtPSAMax`, `AtELossModel`, `AtFitterTask`. + +Default to `-!` unless disk persistence is actually required. + +### C++ vs ROOT Typedefs + +- Non-persisted classes (tasks, algorithms, models): use `bool`, `int`, `double`, `std::string`. +- Persisted data classes only: use `Bool_t`, `Int_t`, `Double_t`, etc. + +### Test Isolation + +Unit tests must not access external files or network resources. Hardcode test data inline. + +## Contributing + +- PRs target the `develop` branch; fast-forward only (no merge commits). +- Commit messages: present imperative mood, ≤72 characters. +- All PRs must pass `clang-format`, `clang-tidy`, and unit tests. diff --git a/AtTools/AtELossBetheBloch.cxx b/AtTools/AtELossBetheBloch.cxx index 91c7e283a..4d36362d4 100644 --- a/AtTools/AtELossBetheBloch.cxx +++ b/AtTools/AtELossBetheBloch.cxx @@ -78,6 +78,24 @@ void AtELossBetheBloch::BuildSpline(double E_min_MeV, double E_max_MeV, int nPoi dXdE[i] = (dedxValues[i] > 0) ? 1.0 / dedxValues[i] : 1e10; fdXdE = tk::spline(energies, dXdE); + + // Build Ω²(E) spline: Bohr range variance accumulated from 0 to E. + // ω²_unit [MeV²/mm] = K · z² · (Z/A) · ρ[g/mm³] · mₑc² + // (kK is in MeV cm²/mol; ×100 → MeV mm²/mol; density g/cm³ → g/mm³ = /1000; combined: /10) + double z = IsElectron() ? 1.0 : fPart_q; + double omega2_unit = kK * z * z * (static_cast(fMat_Z) / fMat_A) * (fDensity / 10.0) * kM_e; + + // Build spline of the integrand dΩ²/dE = 1/|dEdx|³, then integrate segment-by-segment using + // Simpson's rule (exact for cubics, O(h⁴)) rather than the O(h²) trapezoidal rule. + std::vector integrandValues(nPoints); + for (int i = 0; i < nPoints; ++i) + integrandValues[i] = (dedxValues[i] > 0) ? 1.0 / (dedxValues[i] * dedxValues[i] * dedxValues[i]) : 0.0; + tk::spline integrand(energies, integrandValues); + + std::vector rangeVar(nPoints, 0.0); + for (int i = 1; i < nPoints; ++i) + rangeVar[i] = rangeVar[i - 1] + omega2_unit * integrand.integrate(energies[i - 1], energies[i]); + fRangeVariance = tk::spline(energies, rangeVar); } double AtELossBetheBloch::GetdEdx_formula(double energy) const @@ -199,15 +217,12 @@ double AtELossBetheBloch::GetEnergy(double energyIni, double distance) const double AtELossBetheBloch::GetElossStraggling(double energyIni, double energyFin) const { - double dx_mm = GetRange(energyIni, energyFin); - if (dx_mm <= 0) + if (energyIni <= energyFin) 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(fMat_Z) / fMat_A) * fDensity * dx_cm * kM_e; - return std::sqrt(sigma2); + double omega2 = GetRangeVariance(energyIni) - GetRangeVariance(energyFin); + if (omega2 <= 0) + return 0; + return std::abs(GetdEdx(energyFin)) * std::sqrt(omega2); } double AtELossBetheBloch::GetdEdxStraggling(double energyIni, double energyFin) const @@ -215,7 +230,19 @@ double AtELossBetheBloch::GetdEdxStraggling(double energyIni, double energyFin) double dx_mm = GetRange(energyIni, energyFin); if (dx_mm <= 0) return 0; - return GetElossStraggling(energyIni, energyFin) / dx_mm; + double omega2 = GetRangeVariance(energyIni) - GetRangeVariance(energyFin); + if (omega2 <= 0) + return 0; + return std::abs(GetdEdx(energyFin)) * std::sqrt(omega2) / dx_mm; +} + +double AtELossBetheBloch::GetRangeVariance(double energy) const +{ + if (energy <= fRangeVariance.get_x_min()) + return 0; + if (energy >= fRangeVariance.get_x_max()) + return fRangeVariance(fRangeVariance.get_x_max()); + return fRangeVariance(energy); } } // namespace AtTools diff --git a/AtTools/AtELossBetheBloch.h b/AtTools/AtELossBetheBloch.h index a2bcb515c..8dadd222e 100644 --- a/AtTools/AtELossBetheBloch.h +++ b/AtTools/AtELossBetheBloch.h @@ -28,7 +28,8 @@ class AtELossBetheBloch : public AtELossModel { int fMat_A; // target mass number (g/mol) double fI_MeV; // mean excitation energy in MeV - tk::spline fdXdE; // spline of dx/dE vs. energy; integral cached for O(log n) GetRange + tk::spline fdXdE; // spline of dx/dE vs. energy; integral cached for O(log n) GetRange + tk::spline fRangeVariance; // spline of Ω²(E) [mm²] vs. energy; accumulated Bohr range variance public: /** @@ -57,6 +58,7 @@ class AtELossBetheBloch : public AtELossModel { virtual double GetEnergy(double energyIni, double distance) const override; virtual double GetElossStraggling(double energyIni, double energyFin) const override; virtual double GetdEdxStraggling(double energyIni, double energyFin) const override; + virtual double GetRangeVariance(double energy) const override; private: bool IsElectron() const { return std::abs(fPart_mass - kM_e) < 0.01; } diff --git a/AtTools/AtELossBetheBlochTest.cxx b/AtTools/AtELossBetheBlochTest.cxx index 2dfcaec2e..f65ca93b3 100644 --- a/AtTools/AtELossBetheBlochTest.cxx +++ b/AtTools/AtELossBetheBlochTest.cxx @@ -1,5 +1,9 @@ #include "AtELossBetheBloch.h" +#include "AtELossCATIMA.h" + +#include +#include #include #include @@ -145,15 +149,79 @@ TEST_F(AtELossBetheBlochFixture, SettersRebuildSpline) EXPECT_NE(m.GetdEdx(1.0), dedx_before); } -TEST_F(AtELossBetheBlochFixture, BohrStragglingSanity) +// Proton in Ar gas (same material as PiondEdx test). +class AtELossBetheBlochArFixture : public ::testing::Test { +protected: + AtELossBetheBloch model; + // Ar: Z=18, A=40, ρ=1.65e-3 g/cm³, I=188 eV + AtELossBetheBlochArFixture() : model(1.0, kProtonMass, 18, 40, 1.65e-3, 188.0) {} +}; + +TEST_F(AtELossBetheBlochArFixture, RangeVarianceZeroAtLowEnergy) +{ + EXPECT_DOUBLE_EQ(model.GetRangeVariance(0.0), 0.0); + EXPECT_DOUBLE_EQ(model.GetRangeVariance(1e-6), 0.0); +} + +TEST_F(AtELossBetheBlochArFixture, RangeVarianceMonotonicallyIncreasing) +{ + double prev = 0.0; + for (double E : {1.0, 2.0, 5.0, 10.0, 20.0, 50.0}) { + double rv = model.GetRangeVariance(E); + EXPECT_GT(rv, prev) << "RangeVariance not increasing at E=" << E; + prev = rv; + } +} + +TEST_F(AtELossBetheBlochArFixture, StragglingScalesWithSqrtRange) { - double E0 = 5.0; - double Ef = model.GetEnergy(E0, 100.0); - EXPECT_GT(Ef, 0.0); + // For a factor-4 path length ratio, straggling should scale as sqrt(path) → factor ~2. + // Use two different path lengths by picking final energies via GetEnergy. + double E0 = 20.0; + double Ef1 = model.GetEnergy(E0, 500.0); + double Ef2 = model.GetEnergy(E0, 2000.0); // 4× longer path + ASSERT_GT(Ef1, 0.0); + ASSERT_GT(Ef2, 0.0); + + double sigma1 = model.GetElossStraggling(E0, Ef1); + double sigma2 = model.GetElossStraggling(E0, Ef2); + ASSERT_GT(sigma1, 0.0); + ASSERT_GT(sigma2, 0.0); + + // σ scales roughly as sqrt(Δx), so σ2/σ1 ≈ sqrt(4) = 2. + // Allow 20% tolerance: exact scaling only holds when dEdx is constant over the path. + double ratio = sigma2 / sigma1; + EXPECT_NEAR(ratio, 2.0, 0.20 * 2.0); +} + +TEST_F(AtELossBetheBlochArFixture, dEdxStragglingConsistency) +{ + // GetdEdxStraggling * GetRange must equal GetElossStraggling (up to floating point). + double E0 = 10.0; + double Ef = model.GetEnergy(E0, 1000.0); + ASSERT_GT(Ef, 0.0); + + double sigmaE = model.GetElossStraggling(E0, Ef); + double dx = model.GetRange(E0, Ef); + double sigmaDedx = model.GetdEdxStraggling(E0, Ef); + + ASSERT_GT(sigmaE, 0.0); + ASSERT_GT(dx, 0.0); + EXPECT_NEAR(sigmaDedx * dx, sigmaE, 1e-9 * sigmaE); +} +TEST_F(AtELossBetheBlochArFixture, StragglingBroadRange) +{ + // Physical bound: straggling cannot exceed energy loss (σ(ΔE) < ΔE). + double E0 = 10.0; + double Ef = 5.0; + double dE = E0 - Ef; double sigma = model.GetElossStraggling(E0, Ef); + EXPECT_GT(sigma, 0.0); - EXPECT_LT(sigma, E0); // straggling must be less than total energy + EXPECT_LT(sigma, dE); + // Also check the ratio is sub-50% (physically meaningful bound) + EXPECT_LT(sigma / dE, 0.5); } TEST_F(AtELossBetheBlochFixture, BlochApprox) @@ -168,3 +236,67 @@ TEST_F(AtELossBetheBlochFixture, BlochApprox) // Bloch approximation should agree to within 30% of the explicit I value EXPECT_NEAR(rangeBloch, rangeExact, 0.30 * rangeExact); } + +/** + * Benchmark fixture: compare AtELossBetheBloch against AtELossCATIMA configured for pure Bohr + * straggling (z_eff_type::none = bare charge, default calculation = bohr mode). + * + * System: proton in H₂ at 600 Torr (ρ = 6.5643e-5 g/cm³), identical to AtELossCATIMATestFixture. + * Expected agreement: within 10% (residual difference from Lindhard X × γ² ≈ 1 for protons at 1–10 MeV). + */ +class AtELossBetheBlochVsCATIMAFixture : public ::testing::Test { +protected: + static constexpr double kProtonMassAmu = 1.007825031898; + static constexpr double kDensity = 6.5643e-5; // g/cm³ + + AtTools::AtELossBetheBloch bb; + AtTools::AtELossCATIMA catima; + + AtELossBetheBlochVsCATIMAFixture() + : bb(1.0, kProtonMass, 1, 1, kDensity, 19.2), catima(kDensity, catima::Material(1, 1)) + { + catima.SetProjectile(1, 1, kProtonMassAmu); + + // Configure CATIMA to pure Bohr: bare charge (no effective-charge model), + // default calculation=bohr already disables the ATIMA correction. + catima::Config pureBohrCfg; + pureBohrCfg.z_effective = catima::z_eff_type::none; + catima.SetConfig(pureBohrCfg); + } +}; + +TEST_F(AtELossBetheBlochVsCATIMAFixture, CATIMAComparison_RangeVariance) +{ + // Compare accumulated Bohr range variance Ω²(E) [mm²] at three energies. + // Both models implement the same Bohr formula; agreement within 10% is expected. + for (double E : {1.0, 5.0, 10.0}) { + double bbVal = bb.GetRangeVariance(E); + double catimaVal = catima.GetRangeVariance(E); + ASSERT_GT(bbVal, 0.0) << "BB range variance zero at E=" << E; + ASSERT_GT(catimaVal, 0.0) << "CATIMA range variance zero at E=" << E; + EXPECT_NEAR(bbVal, catimaVal, 0.10 * catimaVal) << "Range variance mismatch at E=" << E << " MeV"; + } +} + +TEST_F(AtELossBetheBlochVsCATIMAFixture, CATIMAComparison_EnergyStraggling) +{ + // Reference energy pairs from AtELossCATIMATestFixture::TestEnergyLossStraggling. + // CATIMA (pure Bohr config) is used as the reference; BB must agree within 10%. + const double mass = kProtonMassAmu; + struct Case { + double eIni, eFin; + }; + const Case cases[] = { + {1.0, 0.75 * mass}, // narrow: ~25% energy loss + {5.0, 3.58164 * mass}, // narrow: ~28% energy loss + {5.0, 1.0}, // wide: 80% energy loss + {10.0, 1.0}, // wide: 90% energy loss, exercises Bragg-peak region + }; + for (auto &c : cases) { + double bbSigma = bb.GetElossStraggling(c.eIni, c.eFin); + double catimaSigma = catima.GetElossStraggling(c.eIni, c.eFin); + ASSERT_GT(bbSigma, 0.0) << "BB straggling zero for E0=" << c.eIni; + ASSERT_GT(catimaSigma, 0.0) << "CATIMA straggling zero for E0=" << c.eIni; + EXPECT_NEAR(bbSigma, catimaSigma, 0.10 * catimaSigma) << "Energy straggling mismatch for E0=" << c.eIni << " MeV"; + } +} diff --git a/AtTools/AtELossCATIMA.cxx b/AtTools/AtELossCATIMA.cxx index 3f85a981d..7920bdd58 100644 --- a/AtTools/AtELossCATIMA.cxx +++ b/AtTools/AtELossCATIMA.cxx @@ -23,7 +23,7 @@ double AtELossCATIMA::GetdEdx(double energy) const return 0; } - catima::Result result = catima::calculate(*fProjectile, *fMaterial, energy / fProjectileMassAmu); + catima::Result result = catima::calculate(*fProjectile, *fMaterial, energy / fProjectileMassAmu, fConfig); double dEdx = result.dEdxi * fDensity; // MeV/cm return dEdx / 10.0; // convert to MeV/mm } @@ -38,13 +38,14 @@ double AtELossCATIMA::GetRange(double energyIni, double energyFin) const if (energyFin == 0) { fProjectile->T = energyIni / fProjectileMassAmu; - return catima::range(*fProjectile, *fMaterial) / fDensity * 10.; + return catima::range(*fProjectile, *fMaterial, fConfig) / fDensity * 10.; } double remainingEnergy{energyIni}; double range{0}; while (remainingEnergy > energyFin) { - catima::Result result = catima::calculate(*fProjectile, *fMaterial, remainingEnergy / fProjectileMassAmu); + catima::Result result = + catima::calculate(*fProjectile, *fMaterial, remainingEnergy / fProjectileMassAmu, fConfig); double dEdx = result.dEdxi * fDensity; double DE = dEdx * fRangeStepSize / 10.; @@ -65,7 +66,8 @@ double AtELossCATIMA::GetEnergy(double energyIni, double distance) const double remainingEnergy{energyIni}; double range{0}; while (range < distance) { - catima::Result result = catima::calculate(*fProjectile, *fMaterial, remainingEnergy / fProjectileMassAmu); + catima::Result result = + catima::calculate(*fProjectile, *fMaterial, remainingEnergy / fProjectileMassAmu, fConfig); double dEdx = result.dEdxi * fDensity; double DE{}; @@ -90,7 +92,7 @@ double AtELossCATIMA::GetRangeVariance(double energy) const return 0; } auto range_var = - catima::range_variance(*fProjectile, energy / fProjectileMassAmu, *fMaterial); // range var in (g/cm^2)^2 + catima::range_variance(*fProjectile, energy / fProjectileMassAmu, *fMaterial, fConfig); // range var in (g/cm^2)^2 LOG(debug) << "Range variance in (g/cm^2)^2: " << range_var << " for energy: " << energy / fProjectileMassAmu << " MeV/u"; range_var /= fDensity * fDensity; // convert to (cm)^2 @@ -110,7 +112,7 @@ double AtELossCATIMA::GetElossStraggling(double energyIni, double energyFin) con return 0; } auto energy_strag = catima::energy_straggling_from_E(*fProjectile, energyIni / fProjectileMassAmu, - energyFin / fProjectileMassAmu, *fMaterial); + energyFin / fProjectileMassAmu, *fMaterial, fConfig); return energy_strag; } double AtELossCATIMA::GetdEdxStraggling(double energyIni, double energyFin) const @@ -142,9 +144,10 @@ AtELossCATIMA::GetBraggCurve(double energy, double rangeStepSize, double totalFr double range{}; while (remainingEnergy / energy > totalFractionELoss) { - catima::Result result = catima::calculate(*fProjectile, *fMaterial, remainingEnergy / fProjectileMassAmu); - double dEdx = result.dEdxi * fDensity; - braggCurve.push_back(std::make_pair(dEdx, range)); + catima::Result result = + catima::calculate(*fProjectile, *fMaterial, remainingEnergy / fProjectileMassAmu, fConfig); + double dEdx = result.dEdxi * fDensity; // MeV/cm + braggCurve.push_back(std::make_pair(dEdx / 10.0, range)); // store in MeV/mm double DE = dEdx * rangeStepSize / 10.; if (DE > remainingEnergy) diff --git a/AtTools/AtELossCATIMA.h b/AtTools/AtELossCATIMA.h index bf63862e7..f7b3fbb63 100644 --- a/AtTools/AtELossCATIMA.h +++ b/AtTools/AtELossCATIMA.h @@ -5,6 +5,7 @@ #include "AtELossModel.h" #include +#include #include #include #include @@ -19,6 +20,7 @@ class AtELossCATIMA : public AtELossModel { double fProjectileMassAmu{-1}; /// Mass of the projectile in amu (atomic mass units). double fRangeStepSize{0.1}; // mm + catima::Config fConfig{catima::default_config}; public: /** @@ -68,6 +70,7 @@ class AtELossCATIMA : public AtELossModel { * @param[in] stepSize The step size used for ranges. It must be input in mm. */ void SetRangeStepSize(double stepSize) { fRangeStepSize = stepSize; } + void SetConfig(catima::Config cfg) { fConfig = cfg; } }; } // namespace AtTools diff --git a/docs/contributing/code-style.md b/docs/contributing/code-style.md new file mode 100644 index 000000000..f62af2648 --- /dev/null +++ b/docs/contributing/code-style.md @@ -0,0 +1,68 @@ +# Code Style + +This page covers naming conventions, formatting, and static-analysis tools. Contributor-critical ROOT and persistence rules live in [guide.md](guide.md). + +## Naming Conventions + +### Libraries and Classes + +- **Libraries**: a folder named `AtFoo` produces shared object `libAtFoo` +- **Classes**: non-namespaced classes use the `At` prefix in CamelCase — `AtBaseEvent`, `AtPSAMax`, `AtFitterTask` +- Namespaced code (e.g., `EventFit::`, `kf::`) omits the `At` prefix + +### Data Members + +- Private and protected members: `f` prefix + capital letter — `fTrackID`, `fDriftVelocity` +- Boolean flags: `k` prefix — `kIsGood`, `kInitialized` +- Local variables and function parameters: no prefix, lowercase camelCase + +### Methods + +Follow ROOT naming conventions: +- Getters: `GetXaxis()` not `GetXAxis()` (avoid back-to-back capitals) +- Setters: `SetThreshold()` +- Boolean queries: `IsValid()`, `HasHits()` + +### ClassDef / ClassImp + +Use `ClassDef` and `ClassImp` only for classes that need ROOT-generated class metadata or streamer support. + +For reflection-only classes, prefer a `-!` LinkDef entry. Use `+` only for classes written to disk. + +When a streamed class changes memory layout, increment the version number in `ClassDef`: + +```cpp +ClassDef(AtMyData, 2); // bump when adding/removing/reordering members +``` + +## Formatting + +The project uses `clang-format-17` with the repo-root `.clang-format`. + +Format one file: + +```bash +clang-format-17 -i path/to/file.cxx +``` + +Format changed files: + +```bash +scripts/formatAll.sh +``` + +## Static Analysis + +The repo also carries `.clang-tidy` configuration and CI checks for: + +- `clang-format` +- `clang-tidy` +- `iwyu` + +Run clang-tidy manually on a file: + +```bash +clang-tidy path/to/file.cxx -- -I... +``` + +For framework-specific contributor rules, use [guide.md](guide.md). diff --git a/docs/contributing/guide.md b/docs/contributing/guide.md new file mode 100644 index 000000000..7713eba04 --- /dev/null +++ b/docs/contributing/guide.md @@ -0,0 +1,95 @@ +# Contributor Guide + +## Persisted vs Non-Persisted Code + +Use plain C++ types in non-persisted classes such as tasks, algorithms, models, samplers, and helpers: + +```cpp +bool flag; +int count; +double value; +std::string name; +``` + +Use ROOT typedefs only in classes written to ROOT files: + +```cpp +Bool_t fFlag; +Int_t fCount; +Double_t fValue; +``` + +## ROOT Dictionaries + +Every class with `ClassDef` must appear in the appropriate `*LinkDef.h`. +Classes listed in `*LinkDef.h` do not automatically need `ClassDef` or `ClassImp`. + +Preferred suffix rules for new or edited code: + +- `ClassName +;` + for classes written to disk +- `ClassName -!;` + for reflection-only classes that should not get streamers + +Default to `-!` unless the class is actually persisted. + +If you add a nested type inside a persisted class and ROOT needs to stream it, register that nested type explicitly as well. +If a persisted class owns streamed pointees through pointers or smart pointers, the pointee types still need the required dictionary coverage. + +Current LinkDefs in the tree still contain legacy `+` entries for some non-persisted classes. Treat the rule above as the preferred direction for new work, not as permission for mechanical repo-wide rewrites. + +## Persisted Data Model Conventions + +Do not assume older ROOT folklore applies everywhere in this tree. Persisted classes here already use STL containers and `std::unique_ptr` where ROOT dictionary support is in place. + +Examples in the current branch: + +- `AtEvent` stores hits in `std::vector>` +- `AtTrack` stores hits and pattern state with smart pointers +- `AtTrackingEvent` stores fitted tracks in `std::vector>` + +For new code, follow the local data-model pattern of the owning module. Do not "simplify" these classes back to raw pointers or `TClonesArray` just because older ROOT examples elsewhere do that. + +## Include and Header Norms + +- Public headers are copied into `build/include`, so code in this tree normally includes project headers by basename such as `"AtEvent.h"`, not module-qualified paths. +- Prefer the include style already used in the surrounding module. Avoid mechanical rewrites to `AtData/AtEvent.h`-style paths. +- Public headers commonly use forward declarations for ROOT and FairRoot types. Follow the local header pattern instead of eagerly adding heavyweight includes. + +## When You Add Code + +For a new persisted class, usually update all of: + +- header/source in the owning module +- module `CMakeLists.txt` +- module `*LinkDef.h` +- tests +- docs if branch flow or the runtime data model changes + +For a new non-persisted task or algorithm, usually update all of: + +- header/source in the owning module +- module `CMakeLists.txt` +- module `*LinkDef.h` with `-!` +- tests +- [branch-io-contracts.md](../reference/branch-io-contracts.md) if it reads or writes FairRoot branches + +## CMake and Test Registration + +- `generate_target_and_root_library(...)` + wires libraries and dictionaries together +- `attpcroot_generate_tests(...)` + wires test binaries into CTest + +If you update code but miss one of those registrations, the branch will often build partially and fail later. + +For test writing rules and registration patterns see [testing.md](../tooling/testing.md). + +## ROOT / FairRoot Gotchas + +- Objects read from `FairRootManager` or pulled out of framework-managed `TClonesArray` containers are framework-owned. +- Do not manually `delete` branch objects returned by the framework. +- Do not wrap FairRoot-owned branch objects in smart pointers. +- Macros still use raw `new` for tasks and generators because ownership is transferred through framework registration such as `AddTask(...)` and generator setup. +- Inside tasks and helper classes, prefer the ownership pattern already used locally. In current code that often means `std::unique_ptr` for owned algorithms or models. +- If a ROOT macro cannot resolve a symbol, first suspect the library, dictionary, or environment setup before adding project header includes. diff --git a/docs/contributing/new-module.md b/docs/contributing/new-module.md new file mode 100644 index 000000000..c5c263ad4 --- /dev/null +++ b/docs/contributing/new-module.md @@ -0,0 +1,71 @@ +# Adding a New Module Library + +Each ATTPCROOT module is a separate CMake shared library with an associated ROOT dictionary. + +## CMakeLists.txt Pattern + +```cmake +set(LIBRARY_NAME MyModule) +set(SRCS + MyClass.cxx + AnotherClass.cxx +) +set(DEPENDENCIES + ATTPCROOT::AtData + ROOT::Core + FairRoot::Base +) + +# Optional: unit tests +set(TEST_SRCS MyClassTest.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} +) +``` + +## LinkDef File + +Every class with a `ClassDef` macro in its header must appear in `MyModuleLinkDef.h`. ROOT uses this file to generate the dictionary. + +```cpp +#ifdef __CINT__ +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MyDataClass+; // preferred for persisted classes +#pragma link C++ class MyTask-!; // preferred for reflection-only code +#pragma link C++ class MyAlgorithm-!; + +#endif +``` + +For the preferred streamer suffix rules, persisted-vs-non-persisted type rules, and notes on legacy LinkDef entries already in the tree, use [guide.md](guide.md). + +## Common FairTask Pattern + +Most framework tasks follow: + +- `Init()` + load/register branches and initialize owned algorithms +- `SetParContainers()` + load runtime-db containers when needed +- `Exec()` + clear outputs, read the input `TClonesArray` and current event object, and write the next branch container + +If a task reads or writes documented branches, update [branch-io-contracts.md](../reference/branch-io-contracts.md). + +## Registering with the Build + +Add the new subdirectory to the top-level `CMakeLists.txt`: + +```cmake +add_subdirectory(MyModule) +``` + +And add `ATTPCROOT::MyModule` to the `DEPENDENCIES` of any module that needs it. diff --git a/docs/development/fitting-status.md b/docs/development/fitting-status.md new file mode 100644 index 000000000..bdcdbd78e --- /dev/null +++ b/docs/development/fitting-status.md @@ -0,0 +1,19 @@ +# Fitting Status + +Fitting is unstable on this branch. This page records only current status. + +## Current State + +The current branch contains multiple fitting approaches, including: + +- GenFit-based fitting through `AtFitterTask` and `AtFITTER::AtGenfit` +- Monte-Carlo-based fitting through `AtMCFitterTask` +- Bragg-curve-related fitting logic in other parts of the codebase + +No stable long-term fitting architecture is documented here. + +## Guidance + +- Treat the current branch code as the source of truth. +- Do not assume the current fitters represent the intended long-term design. +- Inspect source directly before editing fitter code. diff --git a/docs/index.md b/docs/index.md new file mode 100644 index 000000000..e8f9c43ab --- /dev/null +++ b/docs/index.md @@ -0,0 +1,40 @@ +# Development Documentation + +Agent-facing docs for the current branch state of ATTPCROOT. + +## Tasks + +- install or first-time configure: [tooling/installation.md](tooling/installation.md) +- set up a session / build / run tests: [tooling/daily-use.md](tooling/daily-use.md) +- write and register tests: [tooling/testing.md](tooling/testing.md) +- add code safely: [contributing/guide.md](contributing/guide.md), [contributing/new-module.md](contributing/new-module.md) +- check formatting rules: [contributing/code-style.md](contributing/code-style.md) +- trace branch flow: [reference/branch-io-contracts.md](reference/branch-io-contracts.md), [subsystems/reconstruction-pipeline.md](subsystems/reconstruction-pipeline.md), [subsystems/simulation-pipeline.md](subsystems/simulation-pipeline.md) +- understand runtime objects: [reference/data-model.md](reference/data-model.md) +- find subsystem ownership: [reference/modules.md](reference/modules.md) +- find an example macro: [reference/macro-cookbook.md](reference/macro-cookbook.md) +- configure detector parameters: [reference/parameters.md](reference/parameters.md) +- prepare detector geometry: [subsystems/geometry.md](subsystems/geometry.md) +- browse data interactively: [subsystems/visualization.md](subsystems/visualization.md) + +## Architecture + +- module map: [reference/modules.md](reference/modules.md) +- data model: [reference/data-model.md](reference/data-model.md) +- simulation flow: [subsystems/simulation-pipeline.md](subsystems/simulation-pipeline.md) +- reconstruction flow: [subsystems/reconstruction-pipeline.md](subsystems/reconstruction-pipeline.md) +- generators: [subsystems/generators.md](subsystems/generators.md) +- PSA: [subsystems/psa.md](subsystems/psa.md) +- energy loss: [subsystems/energy-loss.md](subsystems/energy-loss.md) +- geometry: [subsystems/geometry.md](subsystems/geometry.md) +- visualization: [subsystems/visualization.md](subsystems/visualization.md) +- parameters: [reference/parameters.md](reference/parameters.md) + +## Branch-Specific Notes + +- fitting status: [fitting-status.md](fitting-status.md) + +## Scope Notes + +- These docs describe the current branch state. +- They are framework-centric and keep experiment-specific detail light. diff --git a/docs/reference/branch-io-contracts.md b/docs/reference/branch-io-contracts.md new file mode 100644 index 000000000..71d1e5d97 --- /dev/null +++ b/docs/reference/branch-io-contracts.md @@ -0,0 +1,40 @@ +# Branch and IO Contracts + +Default branch names, contained types, and persistence settings for major tasks. These come from task constructors and `Init()` implementations in the current source. + +Every branch listed here is a `TClonesArray` holding one event object at slot `0` — not a bare pointer: + +```cpp +auto *arr = dynamic_cast(ioMan->GetObject("AtEventH")); +auto *event = dynamic_cast(arr->At(0)); +``` + +## Reconstruction Tasks + +| Task | Input branch / type | Output branch / type | Persist default | Branch setters | +|------|---------------------|----------------------|-----------------|----------------| +| `AtUnpackTask` | external file via unpacker | `AtRawEvent` → `TClonesArray` | **true** | — | +| `AtFilterTask` | `AtRawEvent` → `TClonesArray` | `AtRawEventFiltered` → `TClonesArray` | false | `SetInputBranch`, `SetOutputBranch` | +| `AtPSAtask` | `AtRawEvent` → `TClonesArray` | `AtEventH` → `TClonesArray` | false | `SetInputBranch`, `SetOutputBranch` | +| `AtDataCleaningTask` | `AtEventH` → `TClonesArray` | `AtEventCleaned` → `TClonesArray` | **true** | `SetInputBranch`, `SetOutputBranch` | +| `AtPRAtask` | `AtEventH` → `TClonesArray` | `AtPatternEvent` → `TClonesArray` | false | `SetInputBranch`, `SetOutputBranch` | +| `AtSampleConsensusTask` | `AtEventH` → `TClonesArray` | `AtPatternEvent` → `TClonesArray` | false | `SetInputBranch`, `SetOutputBranch` | +| `AtFitterTask` | `AtPatternEvent` → `TClonesArray` | `AtTrackingEvent` → `TClonesArray` | false | `SetInputBranch`, `SetOutputBranch` | +| `AtMCFitterTask` | `AtPatternEvent` → `TClonesArray` | `AtMCResult` (**true**), `SimEvent` (false), `SimRawEvent` (false) | per-branch | `SetPatternBranchName`, `SetSaveResult/Event/RawEvent` | + +## Simulation and Digitization Tasks + +| Task | Input branch / type | Output branch / type | Persist default | Branch setters | +|------|---------------------|----------------------|-----------------|----------------| +| `AtClusterizeTask` | `AtTpcPoint` → `TClonesArray` | `AtSimulatedPoint` → `TClonesArray` | false | `SetPersistence` | +| `AtPulseTask` | `AtSimulatedPoint` → `TClonesArray`, optional `AtTpcPoint` → `TClonesArray` | `AtRawEvent` → `TClonesArray` (**true**); re-registers `AtTpcPoint` (false) | true / false | `SetOutputBranch`, `SetPersistence`, `SetPersistenceAtTpcPoint`, `SetSaveMCInfo` | + +## Reading This Page + +- These are defaults, not hard-coded universals. Most tasks expose `SetInputBranch` / `SetOutputBranch`. +- `AtPRAtask` and `AtSampleConsensusTask` write the same output branch name — only one should be active per run. +- `AtMCFitterTask` writes three separate branches; each has its own save toggle. +- `AtFittedTrack` (inside `AtTrackingEvent`) is currently unstable — treat layout as best-effort. +- In the AT-TPC path, `AtTpcPoint` is a branch name whose current container holds `AtMCPoint` objects. +- If you add a task that reads or writes FairRoot branches, add a row above. +- For what the contained objects hold (fields, ownership), see [data-model.md](data-model.md). diff --git a/docs/reference/data-model.md b/docs/reference/data-model.md new file mode 100644 index 000000000..2b5366231 --- /dev/null +++ b/docs/reference/data-model.md @@ -0,0 +1,42 @@ +# Data Model + +Core runtime objects passed between simulation, unpacking, reconstruction, and fitting. + +Most FairRoot branches in this tree are not bare `AtEvent*` or `AtRawEvent*` objects. They are `TClonesArray` branch containers, usually holding one event object at slot `0`. Agents should treat the branch name and the contained event class as related but distinct concepts. + +## Reconstruction Objects + +| Object | Typical branch container | Main producer | Main consumer | Persisted | Key downstream fields | +|--------|--------------------------|---------------|---------------|-----------|-----------------------| +| `AtRawEvent` | `TClonesArray` branch such as `AtRawEvent` / `AtRawEventFiltered` | `AtUnpackTask`, `AtPulseTask` | `AtFilterTask`, `AtPSAtask` | yes | pad traces, aux/FPN pads, good flag, optional MC map | +| `AtEvent` | `TClonesArray` branch such as `AtEventH` / `AtEventCleaned` | `AtPSAtask` | `AtDataCleaningTask`, `AtPRAtask`, `AtSampleConsensusTask` | yes | hit list, event charge, mesh signal | +| `AtPatternEvent` | `TClonesArray` branch such as `AtPatternEvent` / `AtPatternEventModified` | `AtPRAtask`, `AtSampleConsensusTask` | `AtFitterTask`, `AtMCFitterTask` | yes | candidate tracks, noise hits | +| `AtTrack` | contained inside `AtPatternEvent`, not its own top-level branch | pattern-recognition code | fitters, track-level analysis | yes | track ID, hit/cluster-hit collections, geometric estimates, Bragg-curve values | +| `AtTrackingEvent` | `TClonesArray` branch `AtTrackingEvent` | `AtFitterTask` | downstream analysis | yes | fitted tracks, optional copied track array, event vertex fields | +| `AtFittedTrack` | contained inside `AtTrackingEvent`, not its own top-level branch | fitter implementations | downstream analysis | yes | current branch layout only; treat as unstable | +| `AtMCResult` | `TClonesArray` branch `AtMCResult` | `AtMCFitterTask` | MC-fitting analysis | yes | objective value and fitted/sampled parameter summary | + +## Simulation Objects + +| Object | Typical branch container | Main producer | Main consumer | Persisted | Notes | +|--------|--------------------------|---------------|---------------|-----------|-------| +| `AtMCPoint` | MC-point branches such as `AtTpcPoint` | Geant/VMC transport, `AtSimpleSimulation` | `AtClusterizeTask` | yes | primary simulation-to-digitization handoff object; also used for MC bookkeeping in simulated-data reconstruction | +| `AtSimulatedPoint` | `TClonesArray` branch `AtSimulatedPoint` | `AtClusterizeTask` | `AtPulseTask` | yes | intermediate ionization-electron / charge-cluster representation between MC points and pad traces | +| `AtMCTrack` | `TClonesArray` branch `MCTrack` | simulation stack | analysis, truth matching | yes | simulated particle tracks | +| `AtVertexPropagator` | no FairRoot branch | generators/runtime simulation logic | generators and downstream simulation code | no | singleton shared state | + +## Container Pattern + +```text +FairRoot branch name -> TClonesArray container -> runtime object + +AtRawEvent -> TClonesArray -> AtRawEvent +AtEventH -> TClonesArray -> AtEvent +AtPatternEvent -> TClonesArray -> AtPatternEvent +AtTrackingEvent -> TClonesArray -> AtTrackingEvent +MCTrack -> TClonesArray -> AtMCTrack +AtTpcPoint -> TClonesArray -> AtMCPoint +AtSimulatedPoint-> TClonesArray -> AtSimulatedPoint +``` + +For the tasks that produce and consume these objects, see [branch-io-contracts.md](branch-io-contracts.md). For pipeline order, see [simulation-pipeline.md](../subsystems/simulation-pipeline.md) and [reconstruction-pipeline.md](../subsystems/reconstruction-pipeline.md). diff --git a/docs/reference/macro-cookbook.md b/docs/reference/macro-cookbook.md new file mode 100644 index 000000000..6e4938fba --- /dev/null +++ b/docs/reference/macro-cookbook.md @@ -0,0 +1,26 @@ +# Macro Cookbook + +Example macro starting points referenced by these docs. + +## Notes + +- `macro/examples/` contains adaptation starting points. +- `macro/tests/` contains repeatable workflow and regression-style examples. +- Many macros are experiment-specific or user-specific and are intentionally omitted. +- Files under `macro/` are ROOT / Cling scripts, not normal compiled translation units. +- ROOT/core headers such as `TClonesArray.h` or `TTreeReader.h`, and STL headers a macro genuinely uses, are normal. +- Do not treat a macro like standalone C++ or try to fix missing symbols by adding project headers first. + +## Macro Do / Don't + +- Do: copy an existing pattern from `macro/examples/` or `macro/tests/`. +- Don't: start by adding `#include "At*.h"` to compensate for missing dictionaries or libraries, or rewrite a macro as if it were a compiled source file. + +## Starting Points + +| Path | Use | Caveat | +|------|-----|--------| +| `macro/examples/run_sim.C` | minimal simulation structure | simple example, not a full workflow template | +| `macro/examples/rundigi_sim.C` | simulation plus downstream task ordering | parts of the task usage are older and should be checked against current APIs | +| `macro/tests/AT-TPC/run_unpack_attpc.C` | unpacking plus PSA/cleaning flow | experiment-facing example, not a framework contract | +| `macro/tests/AT-TPC/run_sim_attpc.C` | fuller simulation test flow | oriented toward test workflow rather than minimal onboarding | diff --git a/docs/reference/modules.md b/docs/reference/modules.md new file mode 100644 index 000000000..919e81a0f --- /dev/null +++ b/docs/reference/modules.md @@ -0,0 +1,30 @@ +# Module Overview + +ATTPCROOT is organized as CMake library targets layered around data, simulation, unpacking, reconstruction, and analysis. + +## Modules + +`Also inspect` points to the most useful overview docs for that module. + +| Module | Role | Main subareas | Common edit reasons | Also inspect | +|--------|------|---------------|---------------------|--------------| +| `AtData` | persisted event and track containers | event, track, pattern classes | add/modify data objects | [data-model.md](data-model.md), [branch-io-contracts.md](branch-io-contracts.md) | +| `AtReconstruction` | PSA, filtering, pattern recognition, fitting tasks | `AtPulseAnalyzer/`, `AtPatternRecognition/`, `AtPatternModification/`, `AtFitter/`, `AtFilter/` | change branch flow or task behavior | [reconstruction-pipeline.md](../subsystems/reconstruction-pipeline.md), [branch-io-contracts.md](branch-io-contracts.md) | +| `AtUnpack` | experimental data to `AtRawEvent` | unpackers, `GETDecoder2/`, `deprecated/` | add file-format/front-end support | [data-model.md](data-model.md) | +| `AtDigitization` | simulation truth to detector signals | clusterization, pulse generation, trigger/space-charge tasks | change MC-to-signal behavior | [simulation-pipeline.md](../subsystems/simulation-pipeline.md), [branch-io-contracts.md](branch-io-contracts.md) | +| `AtTools` | shared utilities | energy loss, kinematics, hit sampling, cleaning, electronic response | algorithmic utilities used across modules | [energy-loss.md](../subsystems/energy-loss.md) | +| `AtSimulationData` | MC truth and simulation shared state | MC objects, stack, propagator | change truth object layout or simulation state | [simulation-pipeline.md](../subsystems/simulation-pipeline.md), [data-model.md](data-model.md) | +| `AtGenerators` | FairRoot generator implementations | reaction and beam generators | change simulation source behavior | [generators.md](../subsystems/generators.md), [simulation-pipeline.md](../subsystems/simulation-pipeline.md) | +| `AtDetectors` | detector geometry and sensitive detectors | detector-specific subdirs, field/passive geometry | detector implementation changes | [geometry.md](../subsystems/geometry.md) | +| `AtAnalysis` | higher-level analysis code | experiment-facing analysis | analysis task changes | [data-model.md](data-model.md), [reconstruction-pipeline.md](../subsystems/reconstruction-pipeline.md) | +| `AtEventDisplay` | event visualization | `AtTabs/`, `AtSidebar/`, legacy display code | visualization changes | [visualization.md](../subsystems/visualization.md), [data-model.md](data-model.md) | + +## Workflow Assets + +| Directory | Role | +|-----------|------| +| `macro/` | example and test macros | +| `geometry/` | detector geometry assets | +| `parameters/` | runtime parameter files | +| `resources/` | shipped data assets such as tables and cross sections | +| `scripts/` | support scripts and mapping-related assets | diff --git a/docs/reference/parameters.md b/docs/reference/parameters.md new file mode 100644 index 000000000..326db7f70 --- /dev/null +++ b/docs/reference/parameters.md @@ -0,0 +1,93 @@ +# Parameter Files + +ATTPCROOT uses FairRoot's parameter system to manage experiment-specific detector configuration. The main container class is `AtDigiPar`, defined in `AtParameter/AtDigiPar.h`. + +## Loading Parameters in a Macro + +```cpp +FairRuntimeDb *rtdb = run->GetRuntimeDb(); +FairParAsciiFileIo *parIo = new FairParAsciiFileIo(); +parIo->open("parameters/ATTPC.myexperiment.par", "in"); +rtdb->setFirstInput(parIo); +rtdb->getContainer("AtDigiPar"); +``` + +Parameter files live in `parameters/`. Each file typically contains an `[AtDigiPar]` section followed by key-value pairs. Pick the file that matches the experiment or create a new one by copying a current experiment file such as `parameters/ATTPC.e20020.par`. + +## AtDigiPar Parameters + +`AtDigiPar::getParams()` requires the following keys: + +### Electromagnetic Fields + +| Parameter | Type | Unit | Description | +|-----------|------|------|-------------| +| `EField` | Double_t | V/m | Longitudinal drift electric field | +| `BField` | Double_t | T | Longitudinal magnetic field | + +### Detector Geometry + +| Parameter | Type | Unit | Description | +|-----------|------|------|-------------| +| `TBEntrance` | Int_t | time buckets | Beam entrance position at detector entrance | +| `ZPadPlane` | Double_t | mm | Position of the micromegas pad plane (detector length) | + +### Gas Physics + +| Parameter | Type | Unit | Description | +|-----------|------|------|-------------| +| `DriftVelocity` | Double_t | cm/µs | Electron drift velocity | +| `EIonize` | Double_t | eV | Effective ionization energy of the fill gas | +| `Fano` | Double_t | — | Fano factor of the fill gas | +| `CoefL` | Double_t | cm²/µs | Longitudinal diffusion coefficient | +| `CoefT` | Double_t | cm²/µs | Transverse diffusion coefficient | +| `GasPressure` | Double_t | torr | Gas pressure | +| `Density` | Double_t | kg/m³ | Gas density | + +### Electronics + +| Parameter | Type | Unit | Description | +|-----------|------|------|-------------| +| `SamplingRate` | Int_t | MHz | GET electronics sampling frequency | +| `GETGain` | Double_t | fC | Gain from GET electronics | +| `PeakingTime` | Int_t | ns | Electronic response peaking time | +| `Gain` | Double_t | — | Average micromegas amplification factor | + +## Accessing Parameters in Code + +Tasks retrieve `AtDigiPar` from `FairRuntimeDb` in their `Init()`: + +```cpp +auto *rtdb = FairRuntimeDb::instance(); +fPar = dynamic_cast(rtdb->getContainer("AtDigiPar")); +``` + +Then in `Exec()`: + +```cpp +double driftVel = fPar->GetDriftVelocity(); // cm/µs +double bField = fPar->GetBField(); // T +``` + +## Example `AtDigiPar` Block + +```text +[AtDigiPar] +BField:Double_t 2.0 +EField:Double_t 5000 +TBEntrance:Int_t 445 +ZPadPlane:Double_t 1000.0 +EIonize:Double_t 15.603 +Fano:Double_t 0.24 +CoefL:Double_t 0.0003 +CoefT:Double_t 0.0003 +GasPressure:Double_t 600 +Density:Double_t 0.197 +DriftVelocity:Double_t 0.87 +Gain:Double_t 1000 +SamplingRate:Int_t 3 +GETGain:Double_t 1000 +PeakingTime:Int_t 720 +``` + +Copy a complete file from `parameters/` and edit the values for the target experiment. diff --git a/docs/subsystems/energy-loss.md b/docs/subsystems/energy-loss.md new file mode 100644 index 000000000..1f1b419af --- /dev/null +++ b/docs/subsystems/energy-loss.md @@ -0,0 +1,59 @@ +# Energy Loss + +ATTPCROOT provides several energy-loss utilities in `AtTools/`. The modern model interface is `AtTools::AtELossModel`; `AtTools::AtELossManager` is a separate legacy lookup-table helper, not the owner or facade for all models. + +## Main Types + +| Type | Role | +|------|------| +| `AtELossModel` | abstract interface for stopping power, range, energy loss, and residual energy | +| `AtELossCATIMA` | CATIMA-backed implementation | +| `AtELossTable` | table-backed implementation, typically from SRIM-style data | +| `AtELossBetheBloch` | standalone Bethe-Bloch helper for simpler analytic use cases | +| `AtELossManager` | older lookup-table utility; do not treat it as the generic `AtELossModel` entry point | + +## `AtELossModel` Interface + +All `AtELossModel` implementations expose the same core API: + +- `GetdEdx(energy)` +- `GetRange(energyIni, energyFin = 0)` +- `GetEnergyLoss(energyIni, distance)` +- `GetEnergy(energyIni, distance)` + +Units in this layer are MeV and mm unless a class-specific note says otherwise. + +## CATIMA + +`AtELossCATIMA` is the main model-backed implementation. It wraps the CATIMA library and requires both material configuration and projectile configuration. + +Typical setup: + +```cpp +using MaterialComp = std::tuple; + +std::vector material = { + {1, 1, 2}, // example stoichiometry entries + {12, 6, 1}, +}; + +AtTools::AtELossCATIMA model(density_g_cm3, material); +model.SetProjectile(projectileA, projectileZ, projectileMassAmu); + +double dE = model.GetEnergyLoss(energyMeV, distanceMm); +double range = model.GetRange(energyMeV); +double residual = model.GetEnergy(energyMeV, distanceMm); +``` + +CATIMA is fetched by CMake if it is not already available locally. It is the only backend in this tree that currently exposes straggling calculations through the `AtELossModel` interface. + +## Table and Analytic Helpers + +- `AtELossTable` follows the same `AtELossModel` interface, but uses precomputed stopping/range tables. +- `AtELossBetheBloch` is a lighter analytic helper and is not a replacement for the full CATIMA-backed workflow when you need straggling or richer material handling. + +## Legacy Lookup-Table Path + +`AtELossManager` predates the `AtELossModel` hierarchy. It still exists in the tree, but it is a separate lookup-table class with its own API such as `GetEnergyLoss(...)`, `GetFinalEnergy(...)`, and `GetDistance(...)`. + +Use it when you specifically need that legacy behavior; do not document or refactor it as if it were the generic model manager for the rest of `AtTools/`. diff --git a/docs/subsystems/generators.md b/docs/subsystems/generators.md new file mode 100644 index 000000000..289d3eb3a --- /dev/null +++ b/docs/subsystems/generators.md @@ -0,0 +1,36 @@ +# Event Generators + +ATTPCROOT uses FairRoot's generator framework. All reaction generators subclass `AtReactionGenerator`, which itself subclasses `FairGenerator`. + +## AtReactionGenerator + +`AtReactionGenerator` is the abstract base for all physics reaction generators in `AtGenerators/`. + +**What to implement:** Subclasses implement `GenerateReaction(FairPrimaryGenerator*)`. This method is responsible for computing the kinematics of the reaction and adding primaries to the stack. + +**What not to touch:** `ReadEvent()` is declared `final`. It handles the beam/reaction event alternation via `AtVertexPropagator` and calls `GenerateReaction()` internally. Do not override it. + +## AtVertexPropagator + +`AtVertexPropagator` is a singleton that passes vertex position, beam momentum, and kinematics between generators and downstream tasks (e.g., digitization tasks that need beam parameters). + +```cpp +AtVertexPropagator::Instance()->GetVertex(); // read vertex +AtVertexPropagator::Instance()->EndEvent(); // called at end of each event to alternate beam/reaction +``` + +**In unit tests:** Call `AtVertexPropagator::Instance()->ResetForTesting()` in `SetUp()` to start from a clean state between tests. + +## Sequential Decays + +Generators can be chained for multi-step decays. Enable on the upstream generator: + +```cpp +generator->SetSequentialDecay(true); +``` + +The second generator in the chain reads the vertex and momentum from `AtVertexPropagator` rather than sampling a fresh beam. + +## Beam/Reaction Alternation + +By default, `ReadEvent()` alternates between inserting a beam event (no reaction) and a full reaction event. This models the AT-TPC's operation where beam tracks are recorded alongside reaction tracks for calibration. The alternation is controlled via `AtVertexPropagator::EndEvent()`. diff --git a/docs/subsystems/geometry.md b/docs/subsystems/geometry.md new file mode 100644 index 000000000..a52ee6581 --- /dev/null +++ b/docs/subsystems/geometry.md @@ -0,0 +1,132 @@ +# Detector Geometry + +Use this file when the task involves detector geometry assets, simulation geometry loading, or geometry-related navigation. + +## Start Here + +- `AtDetectors/AtTpc/AtTpc.cxx` +- `AtDetectors/AtTpc/AtTpc.h` +- `geometry/` +- `macro/examples/run_sim.C` +- `macro/tests/AT-TPC/run_sim_attpc.C` +- `macro/tests/generateGeometry.sh` + +If the task is about pad centers, pad IDs, or pad-plane layout, also inspect: + +- `AtMap/AtTpcMap.cxx` +- `AtMap/AtSpecMATMap.cxx` + +## Runtime Path + +Simulation macros pass a geometry file to `AtTpc` with `SetGeometryFileName()`: + +```cpp +AtTpc *tpc = new AtTpc("ATTPC", kTRUE); +tpc->SetGeometryFileName("ATTPC_v1.2.root"); +run->AddModule(tpc); +``` + +`AtTpc::ConstructGeometry()` dispatches by extension: + +- `.root` -> `ConstructRootGeometry()` +- `.geo` -> ASCII branch exists but is not the practical path here + +For this tree, `AtTpc` geometry is normally loaded from generated ROOT files. + +## Where Geometry Lives + +`geometry/` contains: + +- generator macros such as `ATTPC_v1_2.C`, `ATTPC_He1bar.C`, `GADGET_II.C` +- generated outputs such as `ATTPC_v1.2.root` and `ATTPC_v1.2_geomanager.root` +- `media.geo` for material definitions +- detector-specific helper assets + +Do not assume geometry outputs are missing. This checkout already contains several generated `.root` geometry files. + +## Generation + +Generate a geometry only if the required output is missing or stale. + +Example: + +```bash +cd geometry +root -l -q ATTPC_v1_2.C +``` + +This writes: + +- `ATTPC_v1.2.root` +- `ATTPC_v1.2_geomanager.root` + +Test helper: + +```bash +macro/tests/generateGeometry.sh +``` + +That script generates only a fixed subset of geometries used by tests. It is not a full geometry build step. + +## Detector Variants + +| Script | Detector | Notes | +|--------|----------|-------| +| `ATTPC_v1_2.C` | Full AT-TPC | Used by most unpacking/display macros | +| `ATTPC_v1_1.C` | Full AT-TPC | Older version; some example sim macros still reference this | +| `ATTPC_Proto_v1_0.C` | Prototype AT-TPC | Shorter drift length (500 mm) | +| `ATTPC_He1bar.C` | AT-TPC, He 1 atm | Used by the AT-TPC integration test | +| `ATTPC_He*.C`, `ATTPC_D*.C`, `ATTPC_H*.C` | Gas/pressure variants | Specific gas fill and pressure combinations | +| `SpecMAT_*.C` | SpecMAT | Different pad geometry; scintillator variants available | +| `GADGET_II.C` | GADGET-II | Cylindrical TPC variant | +| `FRIB_TPC_v*.C` | FRIB TPC | | + +## Macro Anchors + +Use real macros, not isolated snippets, when tracing behavior: + +- `macro/examples/run_sim.C` +- `macro/tests/AT-TPC/run_sim_attpc.C` + +Those show the surrounding `FairRunSim` setup, `run->SetMaterials("media.geo")`, and the other modules commonly registered with the TPC geometry such as `AtCave`, `AtMagnet`, and `AtPipe`. + +## Variant Selection + +There is no single geometry file that is clearly canonical across the tree. + +Different workflows use different files: + +- unpacking and display macros often use `ATTPC_v1.2.root` +- older example simulation macros use `ATTPC_v1.1.root` +- the AT-TPC integration test uses `ATTPC_He1bar.root` + +When editing a workflow, trust the specific macro or test in that workflow over generic prose. + +## Materials + +`run->SetMaterials("media.geo")` supplies material definitions used by geometry generation and simulation setup. + +Common modifications to `geometry/media.geo`: + +- **Gas composition**: update the `AtTPC_Gas` entry — change density, atomic composition, and weight fractions +- **Window materials**: modify the entrance window material entry +- **Add a material**: append a new block following the existing format; the file header documents the required fields (`ncomp`, `aw`, `an`, `dens`, `wm`, `sensflag`, `fldflag`, `fld`, `epsil`) + +If `media.geo` changes, regenerate any geometry outputs that depend on it. + +## Geometry vs Mapping + +Keep these separate: + +- transport geometry: `AtDetectors/AtTpc`, `geometry/*.C`, simulation macros +- pad mapping and display-plane geometry: `AtMap/`, visualization, decoder/mapping setup + +If the task mentions pad geometry or lookup tables, start in `AtMap/`, not `AtDetectors/AtTpc`. + +## Navigation Order + +1. Find the macro, task, or test that names the geometry file. +2. Check whether that `.root` file already exists in `geometry/`. +3. Inspect `AtDetectors/AtTpc/AtTpc.cxx` for the runtime load path. +4. Inspect the matching geometry generator macro in `geometry/` if the asset must change. +5. Branch into `AtMap/` only if the task is really about pad mapping rather than transport geometry. diff --git a/docs/subsystems/psa.md b/docs/subsystems/psa.md new file mode 100644 index 000000000..13af6e209 --- /dev/null +++ b/docs/subsystems/psa.md @@ -0,0 +1,50 @@ +# Pulse Shape Analysis (PSA) + +PSA extracts hit positions and charge information from raw pad traces (`AtRawEvent`) to produce an `AtEvent` containing `AtHit` objects. + +## Class Hierarchy + +`AtPSA` is the abstract base class in `AtReconstruction/AtPulseAnalyzer/`. All concrete PSA implementations inherit from it. + +| Class | Hits per pad | Description | +|-------|-------------|-------------| +| `AtPSAMax` | 1 | Places one hit at the time bucket of peak ADC. Charge stored is the peak ADC value; total trace integral also stored on the hit. Optional sub-TB time correction via charge-weighted average of peak region. | +| `AtPSAFull` | 1 or many | Designed for long traces (e.g. beam tracks). Traces with a signal window <50 TBs get one hit at the peak; traces ≥50 TBs are split into 50-TB windows with one hit per window and charge = average ADC in that window. | +| `AtPSASpectrum` | 1 or many | Uses ROOT `TSpectrum` for peak finding; can find multiple peaks per pad. Supports optional background subtraction and interpolation. | +| `AtPSAHitPerTB` | many | One hit per time bucket above threshold. Produces a dense 3D charge distribution with no peak-finding assumptions. Configurable TB range. | +| `AtPSATBAvg` | many | Divides the trace into windows of N TBs (default 5) and creates one hit per window whose average charge exceeds threshold. Hit Z is at window center. Skips saturated pads when a max threshold is set. Can operate on a pad augment instead of the raw ADC. | +| `AtPSADeconv` | 1 | FFT-deconvolves the electronics response function from the raw trace and applies a Butterworth low-pass filter to recover the ionization current. Reconstructed charge saved as `"Qreco"` augment on the pad. Hit placed at charge-distribution weighted average. Requires a response function (`SetResponse`). | +| `AtPSADeconvFit` | 1 | Extends `AtPSADeconv`; instead of the weighted average, fits the deconvolved charge distribution with a Gaussian using the longitudinal diffusion coefficient for improved Z resolution. | +| `AtPSAIterDeconv` | 1 | Extends `AtPSADeconv` with iterative corrections to the reconstructed current. | +| `AtPSASimple2` | — | **Deprecated.** Use `AtPSASpectrum` or `AtPSAMax` instead. | + +## Running PSA + +PSA is applied by `AtPSAtask`, which wraps any `AtPSA` subclass: + +```cpp +auto psa = std::make_unique(); +psa->SetThreshold(45); // ADC threshold; optional +auto psaTask = new AtPSAtask(std::move(psa)); +psaTask->SetInputBranch("AtRawEventFiltered"); // default: "AtRawEvent" +psaTask->SetPersistence(true); +fRun->AddTask(psaTask); +``` + +## Chaining PSA Methods + +`AtPSAComposite` allows running multiple PSA methods in sequence on the same event, for example to handle different pad types with different algorithms: + +```cpp +auto composite = std::make_unique(); +composite->AddPSA(std::make_unique()); +composite->AddPSA(std::make_unique()); +``` + +## Output + +PSA produces an `AtEvent` stored as a branch in the output TTree. Each `AtHit` in the event carries: +- 3D position (x, y, time-bucket → z via drift velocity) +- Charge (meaning depends on algorithm — peak ADC, window average, etc.) +- Trace integral (sum of all ADC values over the trace) +- Pad number and timestamp metadata diff --git a/docs/subsystems/reconstruction-pipeline.md b/docs/subsystems/reconstruction-pipeline.md new file mode 100644 index 000000000..450286968 --- /dev/null +++ b/docs/subsystems/reconstruction-pipeline.md @@ -0,0 +1,57 @@ +# Reconstruction Pipeline + +The reconstruction pipeline converts raw pad traces into hit clouds, candidate tracks, and fitted track results. Each stage is usually a `FairTask` that reads and writes branches through `FairRootManager`. + +In the current tree, those branches are usually `TClonesArray` containers. Tasks typically fetch the branch container from `FairRootManager`, then read or construct the event object at slot `0`. + +## Shared Early Stages + +```text +AtUnpackTask + │ reads raw GET/GRAW, HDF5, or ROOT data + └─ output branch: AtRawEvent -> TClonesArray[AtRawEvent] + ▼ +AtFilterTask [optional] + │ filters raw traces + └─ output branch: AtRawEventFiltered -> TClonesArray[AtRawEvent] + ▼ +AtPSAtask + │ converts traces into hits + └─ output branch: AtEventH -> TClonesArray[AtEvent] + ▼ +AtDataCleaningTask [optional] + │ removes or adjusts hits before pattern recognition + └─ output branch: AtEventCleaned -> TClonesArray[AtEvent] +``` + +The exact branch names and defaults are documented in [branch-io-contracts.md](../reference/branch-io-contracts.md). + +## Pattern Recognition + +Two supported paths produce `AtPatternEvent → TClonesArray[AtPatternEvent]`. Both populate `fHitArray` on each track and `fNoise` on the event with unassigned hits. They differ in which `AtTrack` fields they fill: + +| | `AtPRAtask` | `AtSampleConsensusTask` | +|-|-------------|-------------------------| +| Algorithm | triplet/hierarchical clustering (TriplClust default; `SetPRAlgorithm(n)`) | sample consensus (`SetPatternType`, `SetEstimator`) | +| `fPattern` | always `AtPatternCircle2D` | configurable: `kLine`, `kCircle2D`, `kRay`, `kY` | +| `fGeoRadius`, `fGeoCenter`, `fGeoTheta`, `fGeoPhi` | **filled** — via internal RANSAC after clustering | **not filled** (remain zero) | +| `fHitClusterArray` | populated | not populated | + +The `fGeo*` fields are the momentum-seed inputs for GenFit-based fitters. `AtSampleConsensusTask` leaves them empty, so it requires additional processing (e.g. `AtPatternModificationTask`) before a GenFit-based fitter can run. + +## Fitting + +After pattern recognition, fitting produces the `AtTrackingEvent` branch container. + +Current fit-related tasks include: + +- `AtFitterTask` +- `AtMCFitterTask` + +Fitting is unstable, so this page intentionally stays shallow. See [fitting-status.md](../development/fitting-status.md). + +## Running Reconstruction + +A typical macro creates a `FairRunAna`, adds tasks in order, and calls `Run(...)`. + +For runtime objects, see [data-model.md](../reference/data-model.md). For example macros, see [macro-cookbook.md](../reference/macro-cookbook.md). diff --git a/docs/subsystems/simulation-pipeline.md b/docs/subsystems/simulation-pipeline.md new file mode 100644 index 000000000..e808e922c --- /dev/null +++ b/docs/subsystems/simulation-pipeline.md @@ -0,0 +1,87 @@ +# Simulation Pipeline + +The simulation pipeline converts a generated reaction into MC truth, simulated detector responses, and finally `AtRawEvent` output that can be passed into reconstruction. + +There are currently two ways to generate the simulation-side `AtMCPoint` input: + +- the full FairRoot/VMC transport path used by most detector simulations +- a lighter-weight model-driven path built around `AtSimpleSimulation` + +## Flow + +The overall simulation flow is shared downstream of `AtMCPoint` generation: + +``` +FairPrimaryGenerator + AtReactionGenerator AtSimpleSimulation + │ │ + ▼ ▼ + Geant4 / VMC transport direct AtMCPoint generation + │ │ + └───────────────┬──────────────────┘ + ▼ + AtMCPoint + ▼ +AtClusterizeTask + │ converts MC-point deposits -> ionization electron clusters + ▼ +AtPulseTask + │ drifts electrons and produces simulated pad traces + └─ output branch: AtRawEvent -> TClonesArray[AtRawEvent] +``` + +The only difference between the two paths is how `AtMCPoint` objects are generated. After that, the downstream digitization flow is the same. + +### AtMCPoint Generation + +- **Primary path:** Geant4/VMC transport through the detector geometry produces `AtMCPoint` objects from the generated particles. +- **Secondary path:** `AtDigitization/AtSimpleSimulation.h` advances particles through the active volume in fixed steps, applies an `AtTools::AtELossModel`, and writes `AtMCPoint` objects directly. + +In this tree, `AtSimpleSimulation` uses straight-line propagation. Future versions may extend this step to non-linear tracks and more general transport models. + +### Shared Downstream Stages + +Once `AtMCPoint` objects exist, the remaining stages are shared: + +- `AtClusterizeTask` converts energy deposits into ionization-electron clusters represented as `AtSimulatedPoint` +- `AtPulseTask` drifts those electrons to the pad plane and produces simulated traces in `AtRawEvent` + +## Key Runtime Objects + +- `AtVertexPropagator` + shared simulation-side state between generators and downstream logic +- `AtMCTrack` + simulated particle tracks +- `AtMCPoint` + MC-point type used by digitization logic; also the concrete hit type written by `AtSimpleSimulation` +- `AtRawEvent` + final simulated raw traces used by reconstruction + +See [data-model.md](../reference/data-model.md) for the object-level view and [branch-io-contracts.md](../reference/branch-io-contracts.md) for task branch names. + +## Event Structure + +The FairRoot/VMC generator path represents each physical beam-induced event as two consecutive FairRoot events: + +- Even-indexed event (0, 2, 4, …): beam phase — the beam particle traverses the detector +- Odd-indexed event (1, 3, 5, …): reaction phase — the reaction products are transported + +Events 0+1 form one complete beam-induced event; events 2+3 form the next, and so on. Code that loops over events or checks event indices must account for this pairing. + +**Track IDs:** Within each FairRoot event the beam particle always has `fTrackID = 0`. Reaction products receive subsequent IDs. + +## Required Pieces + +For the FairRoot/VMC path, a simulation run needs: + +- detector geometry +- a configured `FairPrimaryGenerator` with one or more `AtReactionGenerator` subclasses +- the digitization stages `AtClusterizeTask` and `AtPulseTask` +- the experiment parameter set used by digitization + +For the `AtSimpleSimulation` path, the run replaces VMC transport with: + +- an `AtSimpleSimulation` instance +- one or more configured `AtTools::AtELossModel` instances +- the same downstream digitization stages if `AtRawEvent` output is needed + +See [generators.md](generators.md) for generator behavior and [energy-loss.md](energy-loss.md) for the model layer used by `AtSimpleSimulation`. diff --git a/docs/subsystems/visualization.md b/docs/subsystems/visualization.md new file mode 100644 index 000000000..529773ee4 --- /dev/null +++ b/docs/subsystems/visualization.md @@ -0,0 +1,112 @@ +# Data Visualization + +ATTPCROOT includes an interactive event display built on ROOT's TEve framework. The current API is `AtViewerManager` (in `AtEventDisplay/`). + +## Architecture + +`AtViewerManager` is a singleton GUI/controller. It uses `FairRunAna::Instance()` in `AddTask()` and `Init()`, so the macro must create and configure `FairRunAna` before constructing the viewer. Attach **tabs** (display panels) and optionally **tasks** (FairTask subclasses for online re-analysis). Branch names for the three standard data streams — `AtRawEvent`, `AtEvent`, and `AtPatternEvent` — are selected at runtime from dropdown menus in the GUI sidebar. + +The viewer requires a pad map (`AtMap`) loaded from an experiment-specific XML mapping file. + +## Basic Setup + +```cpp +// Create and configure the FairRoot run before constructing the viewer +FairRunAna *fRun = new FairRunAna(); +fRun->SetSource(new FairFileSource(inputFile)); +fRun->SetSink(new FairRootFileSink(outputFile)); +fRun->SetGeomFile(geoFile); + +// Load the pad mapping +TString mapDir = TString(getenv("VMCWORKDIR")) + "/scripts/e12014_pad_mapping.xml"; +auto fMap = std::make_shared(); +fMap->ParseXMLMap(mapDir.Data()); + +// Create viewer and add tabs +AtViewerManager *eveMan = new AtViewerManager(fMap); + +auto tabMain = std::make_unique(); +tabMain->SetMultiHit(100); +eveMan->AddTab(std::move(tabMain)); + +eveMan->Init(); +``` + +## Adding FairTasks for Online Re-analysis + +`AddTask` registers a `FairTask` that runs each time an event is loaded or re-analyzed (e.g. when PSA parameters are changed interactively): + +```cpp +eveMan->AddTask(new AtPSAtask(std::move(psa))); +``` + +Combined with the PSA sidebar addons (`AtSidebarPSA`, `AtSidebarPSADeconv`, etc.), this allows live adjustment of PSA parameters and immediate re-display without restarting the macro. + +## Available Tabs + +### `AtTabMain` +The primary view. Shows a 3D hit display (TEve), a pad plane heat map, and the waveform of the currently selected pad. Selecting a pad in the pad plane updates the waveform display. + +```cpp +auto tab = std::make_unique(); +tab->SetThreshold(20); // minimum charge to draw a hit +tab->SetMultiHit(100); // max hits in a single pad before that pad is suppressed +tab->SetDrawProjection(true); // draw AtPattern projections on the pad plane +tab->SetHitAttributes(TAttMarker(kPink, kFullDotMedium, 1)); +``` + +### `AtTabBraggCurve` +Extends `AtTabMain`. Adds a Bragg curve panel (dE/dx vs. range) alongside the 3D view, populated from `AtTrack::BraggCurve` data. + +### `AtTabPad` +A configurable grid of pad trace panels. Each cell in the grid can independently display a different data source for the currently selected pad. + +```cpp +// 2-row × 4-column grid +auto tab = std::make_unique(2, 4, "MyTraces"); +tab->DrawADC(0, 0); // pedestal-subtracted ADC in cell (0,0) +tab->DrawRawADC(0, 1); // raw ADC +tab->DrawArrayAug("Qreco", 0, 2); // named array augment (e.g. deconvolved charge) +tab->DrawAuxADC("MeshSignal", 0, 3); // auxiliary pad by name +tab->DrawFPN(1, 0); // fixed pattern noise trace +tab->DrawGenTrace(3, 1, 1); // generic trace by ID +tab->DrawHits(0, 0); // overlay hit markers on the ADC trace +eveMan->AddTab(std::move(tab)); +``` + +### `AtTabMacro` +A canvas-based tab driven entirely by user-supplied lambda functions. Callbacks fire at three granularities: + +```cpp +auto tab = std::make_unique(1, 2, "Custom"); + +// Called once when the tree is opened +tab->SetDrawTreeFunction([](AtTabInfo *info) { /* fill histograms */ }, 0, 0); + +// Called each time an event is loaded +tab->SetDrawEventFunction([](AtTabInfo *info) { /* update per-event plots */ }, 0, 1); + +// Called each time the selected pad changes +tab->SetDrawPadFunction([](AtTabInfo *info, Int_t padNum) { /* pad-level plots */ }, 0, 0); + +eveMan->AddTab(std::move(tab)); +``` + +### `AtTabFission` +Extends `AtTabMain`. Adds hit overlays for fission fragment analysis (uncorrected and corrected hit sets). Exposes `GetFissionBranch()` for branch selection. + +### `AtTabEnergyLoss` +Canvas-based tab that calculates dE/dx for two fission fragments from an `AtFissionEvent`. Stacked histograms, charge-sum and Gaussian-fit methods, and the ratio between the two fragments. + +```cpp +AtTabFission *fissionTab = ...; // add fission tab first +auto elTab = std::make_unique(fissionTab->GetFissionBranch()); +eveMan->AddTab(std::move(elTab)); +``` + +## Notes + +- **Branch selection is in the GUI**, not in the macro. The sidebar dropdowns populate from whatever branches are present in the input file. +- `AtViewerManager::Instance()` returns the singleton after construction. +- Requires a display (X11 or native windowing); cannot run headlessly. +- Pad mapping XML files live in `scripts/`. The correct file is experiment-specific. diff --git a/docs/tooling/daily-use.md b/docs/tooling/daily-use.md new file mode 100644 index 000000000..8746256b8 --- /dev/null +++ b/docs/tooling/daily-use.md @@ -0,0 +1,61 @@ +# Daily Use + +## Before Every Session + +Source the generated environment script from the repo root: + +```bash +source build/config.sh +``` + +This sets the following environment variables: + +| Variable | Purpose | +|----------|---------| +| `ROOTSYS` | ROOT installation prefix | +| `VMCWORKDIR` | Repo root (used by macros to find geometry, parameters, etc.) | +| `ROOT_INCLUDE_PATH` | Headers exposed to the ROOT interpreter | +| `LD_LIBRARY_PATH` | Runtime library search path | +| Geant4 data paths | `G4...DATA` variables needed by Geant4 | + +VSCode terminals auto-source this via the workspace settings. + +If `build/config.sh` does not exist, the project has not been configured yet — see [installation.md](installation.md). + +## Building + +```bash +cmake --build build -j10 +``` + +Adjust `-j` to your CPU count. To rebuild a single module: + +```bash +cmake --build build --target AtSimulationData -j10 +``` + +## Running Tests + +Unit tests are built by default (`BUILD_TESTS=ON`). Binaries land in `build/tests/`. + +Run all tests: + +```bash +cd build && ctest -V +``` + +Run a specific test binary directly: + +```bash +./build/tests/AtSimulationDataTests +./build/tests/AtGeneratorsTests +./build/tests/AtToolsTests +``` + +See [testing.md](testing.md) for how to write and register new tests. + +## Verifying the Environment + +```bash +root-config --version # prints ROOT version +``` diff --git a/docs/tooling/installation.md b/docs/tooling/installation.md new file mode 100644 index 000000000..94c6091f4 --- /dev/null +++ b/docs/tooling/installation.md @@ -0,0 +1,50 @@ +# Installation (First-Time Setup) + +This page covers the one-time steps to configure and build ATTPCROOT on a new machine. After completing these steps, day-to-day work uses only `source build/config.sh` and `cmake --build build` — the environment variables below are never needed again. + +## Prerequisites + +FairSoft and FairRoot must already be installed on the machine. FairSoft provides ROOT, Geant4, and VMC. FairRoot provides the task/run framework. + +## Step 1 — Set Configure-Time Variables + +These two variables tell CMake where FairRoot and FairSoft are installed. They only need to be set in the shell session where you run `cmake` for the first time: + +```bash +export FAIRROOTPATH=/path/to/your/fairroot/install +export SIMPATH=/path/to/your/fairsoft/install +``` + +These are machine-specific. Ask your sysadmin or check the FairRoot/FairSoft install documentation for the correct paths on your system. + +## Step 2 — Configure + +Run CMake out-of-source from the repo root: + +```bash +cmake -S . -B build +``` + +If HDF5 is not bundled with your FairSoft installation, add its prefix: + +```bash +cmake -S . -B build -DCMAKE_PREFIX_PATH=/path/to/hdf5 +``` + +A successful configure generates `build/config.sh`. This file captures the full runtime environment and is sourced before any subsequent build or analysis session. + +## Step 3 — Build + +```bash +cmake --build build -j10 +``` + +Adjust `-j` to your CPU count. + +## After Installation + +From this point on: + +- **Before any session and to rebuild:** see [daily-use.md](daily-use.md) + +The `FAIRROOTPATH` and `SIMPATH` variables from Step 1 are baked into `build/config.sh` and do not need to be set again. diff --git a/docs/tooling/testing.md b/docs/tooling/testing.md new file mode 100644 index 000000000..fefe3264d --- /dev/null +++ b/docs/tooling/testing.md @@ -0,0 +1,80 @@ +# Testing + +ATTPCROOT has two categories of tests: unit tests (compiled, fast) and integration tests (interpreted ROOT macros). + +## Unit Tests + +Unit tests use [Google Test](https://github.com/google/googletest). They live alongside the source files they test — for example, `AtSimulationData/AtVertexPropagatorTest.cxx` tests `AtVertexPropagator`. + +### Writing a Test + +Create a `*Test.cxx` file in the module directory and register it in the module's `CMakeLists.txt`: + +```cmake +set(TEST_SRCS MyClassTest.cxx) +attpcroot_generate_tests(${LIBRARY_NAME}Tests SRCS ${TEST_SRCS} DEPS ${LIBRARY_NAME}) +``` + +`attpcroot_generate_tests` links against Google Test and registers the binary with CTest. The test binary is placed in `build/tests/`. + +Typical binary names follow the module, for example: + +- `AtDataTests` +- `AtGeneratorsTests` +- `AtToolsTests` +- `OpenKFTests` for Eigen/OpenKF-specific reconstruction tests + +### Rules + +- **Tests must not access external files or network resources.** Hardcode small data values inline or use files that ship with the repo under `$VMCWORKDIR`. Tests that require external state are fragile and break on other machines. +- Tests should be fast. Unit tests that take more than a second are a smell. +- For naming conventions see the [wiki](https://github.com/ATTPC/ATTPCROOTv2/wiki/Coding-conventions#test-names). + +### Running Unit Tests + +```bash +cd build && ctest -V # all tests +./build/tests/AtToolsTests # one binary directly +``` + +## Integration Tests + +More involved tests (those that require data files or test full pipelines) live in `macro/tests/`. They are written as ROOT macros (`.C` files) and run interpreted, not compiled. + +These are listed and run separately from unit tests. Use them for: + +- end-to-end workflow checks +- simulation plus reconstruction flows +- experiment-facing regression scenarios + +### Prerequisites + +Integration tests require geometry ROOT files that are **not** committed to the repo. Generate them first: + +```bash +macro/tests/generateGeometry.sh +``` + +Some tests also require external data files that are too large for the repo: + +- **GADGETII tests**: require the fishtank dataset (available on FRIB/NSCL cluster only) +- **SpecMAT tests**: require `TTreesGETrun_9993.root` placed in `macro/tests/SpecMAT/data/` + +### Running All Integration Tests + +```bash +macro/tests/runAllTest.sh +``` + +Results are logged in each test subdirectory as `test.log`. + +See [macro-cookbook.md](macro-cookbook.md) for curated macro starting points. + +## Continuous Integration + +Every PR runs: +1. `clang-format` check +2. `clang-tidy` and `iwyu` static analysis +3. Full unit test build and run (`ctest`) + +Run the tests locally before submitting a PR to catch failures early. diff --git a/macro/tests/ELossComparison/ELossComparison.C b/macro/tests/ELossComparison/ELossComparison.C new file mode 100644 index 000000000..43b19be63 --- /dev/null +++ b/macro/tests/ELossComparison/ELossComparison.C @@ -0,0 +1,477 @@ +/** + * ELossComparison.C + * + * Comparative plots of the three energy loss models in ATTPCROOTv2: + * - AtELossBetheBloch (analytic PDG Bethe-Bloch formula) + * - AtELossCATIMA (CATIMA library wrapper) + * - AtELossTable (SRIM lookup table with spline interpolation) + * + * Four scenarios are plotted: + * 1. Proton in H2 gas @ 600 Torr (all three models) + * 2. Deuteron in D2 gas @ 600 Torr (BB + CATIMA + Table) + * 3. Alpha (4He) in D2 gas @ 600 Torr (BB + CATIMA + Table) + * 4. Proton in He gas @ 700 Torr (BB + CATIMA + Table) + * + * Each scenario produces a canvas with four pads: + * Top-left dE/dx vs kinetic energy (log-log) + * Top-right Range vs kinetic energy (log-log) + * Bottom-left Bragg curve (dE/dx vs distance from entry) + * Bottom-right Range straggling sigma vs kinetic energy (log-log) + * Note: Table straggling requires a full SRIM file with the + * "Multiply" conversion section. Simple SRIM tables (no + * conversion header) are loaded without straggling data and + * are omitted from the straggling pad. + * + * Usage (after sourcing build/config.sh): + * root -l -b -q ELossComparison.C + * + * Output: ELossComparison.pdf (four pages, one per scenario) + * + * Model notes + * ----------- + * "Bethe-Bloch (Bohr)" — AtELossBetheBloch: pure PDG 2022 Eq. 34.1, bare projectile charge, + * Bohr range-straggling integral (dω²/dE = K·z²·(Z/A)·ρ·mₑ / |dEdx|³), √E extrapolation + * below the BB validity threshold. + * + * "CATIMA (default)" — AtELossCATIMA with catima::default_config: + * - Effective charge: Pierce-Blann (z_eff → 0 for protons at low β — suppresses dEdx and + * straggling proportionally to z_eff²). + * - dEdx: Bethe-Bloch + shell corrections + Barkas term + density effect + Lindhard correction; + * SRIM-85 tables at low energies. + * - Straggling: Bohr × Lindhard-X factor (quantum correction, always active); Firsov formula + * used as upper cap below 30 MeV/u — this is LOWER than pure Bohr and causes the straggling + * to drop steeply toward zero at low energy. + * + * "CATIMA (bare, no corr)" — AtELossCATIMA with custom config: + * - Effective charge: none (bare Z, matching BB). + * - dEdx corrections disabled: no Barkas, no Lindhard, no shell corrections. + * - Straggling: still Bohr × Lindhard-X + Firsov cap (not removable via config in the bundled + * CATIMA version). Low-energy dEdx still from SRIM-85 tables, not √E. + * This is the closest CATIMA can get to a pure BB+Bohr calculation; any remaining gap versus the + * red BB curve is due to the Firsov straggling correction and SRIM-85 vs √E extrapolation. + */ + +// --------------------------------------------------------------------------- +// Unit conversion helpers (energy -> MeV, length -> mm) +// --------------------------------------------------------------------------- +double UnitToMeV(const std::string &u) +{ + if (u == "eV") + return 1e-6; + if (u == "keV") + return 1e-3; + if (u == "MeV") + return 1.0; + if (u == "GeV") + return 1e3; + return 0.0; +} + +double UnitToMM(const std::string &u) +{ + if (u == "um") + return 1e-3; + if (u == "mm") + return 1.0; + if (u == "cm") + return 10.0; + if (u == "m") + return 1000.0; + return 0.0; +} + +// --------------------------------------------------------------------------- +// LoadSimpleSrimTable +// +// Load a SRIM-format stopping table that has the standard column layout but +// lacks the "Target Density" header and "Multiply" conversion section. +// +// SRIM stores dE/dx in MeV/(mg/cm2). With known density rho (g/cm3): +// dEdx [MeV/mm] = dEdx [MeV/(mg/cm2)] x rho [g/cm3] x 100 +// +// Range straggling cannot be loaded because AtELossTable::LoadRangeVariance +// is private; the returned model has no straggling data. +// --------------------------------------------------------------------------- +AtTools::AtELossTable *LoadSimpleSrimTable(const std::string &fileName, double density) +{ + std::ifstream file(fileName); + if (!file.is_open()) { + ::Error("LoadSimpleSrimTable", "Cannot open file: %s", fileName.c_str()); + return nullptr; + } + + const double conversion = density * 100.0; // MeV/(mg/cm2) -> MeV/mm + + std::vector energy; + std::vector dEdX; + + std::string line; + while (std::getline(file, line)) { + std::istringstream iss(line); + std::vector tok; + std::string t; + while (iss >> t) + tok.push_back(t); + + // Need: energy Eunit dEdx_e dEdx_n range Runit strag Sunit (>=8 tokens) + if (tok.size() < 8) + continue; + + double en, dedxE, dedxN; + try { + en = std::stod(tok[0]); + dedxE = std::stod(tok[2]); + dedxN = std::stod(tok[3]); + } catch (...) { + continue; + } + + double enMeV = en * UnitToMeV(tok[1]); + double dedxMeVmm = (dedxE + dedxN) * conversion; + + if (enMeV <= 0.0 || dedxMeVmm <= 0.0) + continue; + + energy.push_back(enMeV); + dEdX.push_back(dedxMeVmm); + } + + if (energy.empty()) { + ::Error("LoadSimpleSrimTable", "No valid rows parsed from %s", fileName.c_str()); + return nullptr; + } + + return new AtTools::AtELossTable(energy, dEdX, density); +} + +// --------------------------------------------------------------------------- +// Graph-filling helpers: sample a model over a log-spaced energy grid +// --------------------------------------------------------------------------- +TGraph *MakeDedxGraph(AtTools::AtELossModel *m, double eMin, double eMax, int nPts = 200) +{ + auto *g = new TGraph(nPts); + for (int i = 0; i < nPts; ++i) { + double e = eMin * TMath::Power(eMax / eMin, double(i) / (nPts - 1)); + g->SetPoint(i, e, m->GetdEdx(e)); + } + return g; +} + +TGraph *MakeRangeGraph(AtTools::AtELossModel *m, double eMin, double eMax, int nPts = 200) +{ + auto *g = new TGraph(nPts); + for (int i = 0; i < nPts; ++i) { + double e = eMin * TMath::Power(eMax / eMin, double(i) / (nPts - 1)); + g->SetPoint(i, e, m->GetRange(e)); + } + return g; +} + +// hasStrag: pass false for table models loaded without straggling data +TGraph *MakeStragGraph(AtTools::AtELossModel *m, double eMin, double eMax, bool hasStrag, int nPts = 200) +{ + auto *g = new TGraph(); + if (!hasStrag) + return g; + for (int i = 0; i < nPts; ++i) { + double e = eMin * TMath::Power(eMax / eMin, double(i) / (nPts - 1)); + double s = m->GetRangeStraggling(e); + if (s > 0.0) + g->AddPoint(e, s); + } + return g; +} + +TGraph *MakeBraggGraph(AtTools::AtELossModel *m, double energy, double stepMM = 0.5) +{ + auto pts = m->GetBraggCurve(energy, stepMM); + auto *g = new TGraph(int(pts.size())); + for (int i = 0; i < int(pts.size()); ++i) + g->SetPoint(i, pts[i].second, pts[i].first); + return g; +} + +// --------------------------------------------------------------------------- +// ModelEntry: bundles a model pointer with display properties +// --------------------------------------------------------------------------- +struct ModelEntry { + AtTools::AtELossModel *model; + std::string name; + int color; + int lineStyle; + bool hasStraggling; +}; + +// --------------------------------------------------------------------------- +// DrawScenario: fills a pre-divided 4-pad canvas for one particle/material +// --------------------------------------------------------------------------- +void DrawScenario(const std::string &title, const std::vector &entries, double eMin, double eMax, + double braggEnergy, TCanvas *c) +{ + c->SetTitle(title.c_str()); + c->Divide(2, 2, 0.005, 0.005); + + auto StyleGraph = [](TGraph *g, int col, int sty) { + g->SetLineColor(col); + g->SetLineStyle(sty); + g->SetLineWidth(2); + }; + + // --- Pad 1: dE/dx vs E --- + { + c->cd(1)->SetLogx(); + c->cd(1)->SetLogy(); + c->cd(1)->SetGrid(); + auto *mg = new TMultiGraph(); + auto *leg = new TLegend(0.55, 0.65, 0.88, 0.88); + leg->SetBorderSize(0); + for (auto &e : entries) { + auto *g = MakeDedxGraph(e.model, eMin, eMax); + StyleGraph(g, e.color, e.lineStyle); + mg->Add(g, "L"); + leg->AddEntry(g, e.name.c_str(), "L"); + } + mg->SetTitle(Form("%s;Kinetic energy (MeV);dE/dx (MeV/mm)", title.c_str())); + mg->Draw("A"); + leg->Draw(); + } + + // --- Pad 2: Range vs E --- + { + c->cd(2)->SetLogx(); + c->cd(2)->SetLogy(); + c->cd(2)->SetGrid(); + auto *mg = new TMultiGraph(); + auto *leg = new TLegend(0.15, 0.65, 0.48, 0.88); + leg->SetBorderSize(0); + for (auto &e : entries) { + auto *g = MakeRangeGraph(e.model, eMin, eMax); + StyleGraph(g, e.color, e.lineStyle); + mg->Add(g, "L"); + leg->AddEntry(g, e.name.c_str(), "L"); + } + mg->SetTitle(";Kinetic energy (MeV);Range (mm)"); + mg->Draw("A"); + leg->Draw(); + } + + // --- Pad 3: Bragg curve --- + { + c->cd(3)->SetGrid(); + auto *mg = new TMultiGraph(); + auto *leg = new TLegend(0.15, 0.65, 0.48, 0.88); + leg->SetBorderSize(0); + for (auto &e : entries) { + auto *g = MakeBraggGraph(e.model, braggEnergy); + StyleGraph(g, e.color, e.lineStyle); + mg->Add(g, "L"); + leg->AddEntry(g, e.name.c_str(), "L"); + } + mg->SetTitle(Form("Bragg curve (E_{0} = %.1f MeV);Distance from entry (mm);dE/dx (MeV/mm)", braggEnergy)); + mg->Draw("A"); + leg->Draw(); + } + + // --- Pad 4: Range straggling vs E --- + { + c->cd(4)->SetLogx(); + c->cd(4)->SetLogy(); + c->cd(4)->SetGrid(); + auto *mg = new TMultiGraph(); + auto *leg = new TLegend(0.15, 0.65, 0.48, 0.88); + leg->SetBorderSize(0); + bool anyDrawn = false; + for (auto &e : entries) { + auto *g = MakeStragGraph(e.model, eMin, eMax, e.hasStraggling); + if (g->GetN() == 0) { + delete g; + continue; + } + StyleGraph(g, e.color, e.lineStyle); + mg->Add(g, "L"); + leg->AddEntry(g, e.name.c_str(), "L"); + anyDrawn = true; + } + if (anyDrawn) { + mg->SetTitle(";Kinetic energy (MeV);Range straggling #sigma (mm)"); + mg->Draw("A"); + } + leg->Draw(); + } +} + +// =========================================================================== +// Main macro +// =========================================================================== +void ELossComparison() +{ + TString vmcDir = getenv("VMCWORKDIR") ? getenv("VMCWORKDIR") : "../../.."; + TString resDir = vmcDir + "/resources/energy_loss"; + + // ------------------------------------------------------------------------- + // Physical parameters + // ------------------------------------------------------------------------- + + // Particle rest masses (MeV/c2) for Bethe-Bloch + const double kMproton = 938.272; + const double kMdeuteron = 1875.61; + const double kMalpha = 3727.38; + + // Particle masses in amu for CATIMA + const double kMproton_amu = 1.00728; + const double kMdeuteron_amu = 2.01355; + const double kMalpha_amu = 4.00151; + + // Gas densities (g/cm3) at stated pressures, room temperature (~20 C) + const double kRho_H2_600 = 6.5643e-5; // H2 @ 600 Torr (CATIMA/LISE reference tests) + const double kRho_D2_600 = 1.313e-4; // D2 @ 600 Torr (scaled by molar-mass ratio from H2) + const double kRho_He_700 = 1.53e-4; // He @ 700 Torr (He STP x 700/760 x 273/293) + + // Mean excitation energies (eV) -- PDG 2022 Table 34.1 + const double kI_H = 19.2; + const double kI_He = 41.8; + + // ========================================================================= + // Scenario 1: Proton in H2 @ 600 Torr + // All three models. HinH.txt has the full SRIM format (density + Multiply) + // so the table model includes straggling data. + // ========================================================================= + auto bb_p_H2 = new AtTools::AtELossBetheBloch(1.0, kMproton, 1, 1, kRho_H2_600, kI_H); + + auto catima_p_H2 = + new AtTools::AtELossCATIMA(kRho_H2_600, std::vector>{{1, 1, 1}}); + catima_p_H2->SetProjectile(1, 1, kMproton_amu); + + auto table_p_H2 = new AtTools::AtELossTable(0); + table_p_H2->LoadSrimTable(Form("%s/HinH.txt", resDir.Data())); + table_p_H2->SetDensity(kRho_H2_600); + + // ========================================================================= + // Scenario 2: Deuteron in D2 @ 600 Torr + // ========================================================================= + auto bb_d_D2 = new AtTools::AtELossBetheBloch(1.0, kMdeuteron, 1, 2, kRho_D2_600, kI_H); + + auto catima_d_D2 = + new AtTools::AtELossCATIMA(kRho_D2_600, std::vector>{{2, 1, 1}}); + catima_d_D2->SetProjectile(2, 1, kMdeuteron_amu); + + auto table_d_D2 = LoadSimpleSrimTable(Form("%s/deuteron_D2_600torr.txt", resDir.Data()), kRho_D2_600); + + // ========================================================================= + // Scenario 3: Alpha (4He) in D2 @ 600 Torr + // ========================================================================= + auto bb_a_D2 = new AtTools::AtELossBetheBloch(2.0, kMalpha, 1, 2, kRho_D2_600, kI_H); + + auto catima_a_D2 = + new AtTools::AtELossCATIMA(kRho_D2_600, std::vector>{{2, 1, 1}}); + catima_a_D2->SetProjectile(4, 2, kMalpha_amu); + + auto table_a_D2 = LoadSimpleSrimTable(Form("%s/alpha_D2_600torr.txt", resDir.Data()), kRho_D2_600); + + // ========================================================================= + // Scenario 4: Proton in He @ 700 Torr + // ========================================================================= + auto bb_p_He = new AtTools::AtELossBetheBloch(1.0, kMproton, 2, 4, kRho_He_700, kI_He); + + auto catima_p_He = + new AtTools::AtELossCATIMA(kRho_He_700, std::vector>{{4, 2, 1}}); + catima_p_He->SetProjectile(1, 1, kMproton_amu); + + auto table_p_He = LoadSimpleSrimTable(Form("%s/proton_He_700torr.txt", resDir.Data()), kRho_He_700); + + // ========================================================================= + // "CATIMA bare" config: bare charge + no dEdx corrections — closest analog + // to the pure Bethe-Bloch model available within the CATIMA library. + // Note: the Lindhard-X / Firsov straggling correction is always active in + // the bundled CATIMA source and cannot be disabled via Config. + // ========================================================================= + catima::Config pureCfg; + pureCfg.z_effective = catima::z_eff_type::none; // bare Z, not Pierce-Blann + pureCfg.corrections = + catima::no_barkas | catima::no_lindhard | catima::no_shell_correction; + + auto catimaPure_p_H2 = + new AtTools::AtELossCATIMA(kRho_H2_600, std::vector>{{1, 1, 1}}); + catimaPure_p_H2->SetProjectile(1, 1, kMproton_amu); + catimaPure_p_H2->SetConfig(pureCfg); + + auto catimaPure_d_D2 = + new AtTools::AtELossCATIMA(kRho_D2_600, std::vector>{{2, 1, 1}}); + catimaPure_d_D2->SetProjectile(2, 1, kMdeuteron_amu); + catimaPure_d_D2->SetConfig(pureCfg); + + auto catimaPure_a_D2 = + new AtTools::AtELossCATIMA(kRho_D2_600, std::vector>{{2, 1, 1}}); + catimaPure_a_D2->SetProjectile(4, 2, kMalpha_amu); + catimaPure_a_D2->SetConfig(pureCfg); + + auto catimaPure_p_He = + new AtTools::AtELossCATIMA(kRho_He_700, std::vector>{{4, 2, 1}}); + catimaPure_p_He->SetProjectile(1, 1, kMproton_amu); + catimaPure_p_He->SetConfig(pureCfg); + + // ========================================================================= + // Color / line-style scheme + // ========================================================================= + const int kColBB = kRed + 1; + const int kColCATIMA = kBlue + 1; + const int kColTable = kGreen + 2; + const int kColCATIMApure = kMagenta + 1; + + // ========================================================================= + // Draw and save + // ========================================================================= + auto *c1 = new TCanvas("c1", "Proton in H2", 1200, 900); + { + std::vector ents = { + {bb_p_H2, "Bethe-Bloch (Bohr)", kColBB, 1, true}, + {catima_p_H2, "CATIMA (default)", kColCATIMA, 2, true}, + {table_p_H2, "SRIM Table", kColTable, 3, true}, + {catimaPure_p_H2, "CATIMA (bare, no corr)", kColCATIMApure, 7, true}, + }; + DrawScenario("Proton in H_{2} (600 Torr)", ents, 0.1, 10.0, 5.0, c1); + } + + auto *c2 = new TCanvas("c2", "Deuteron in D2", 1200, 900); + { + std::vector ents = { + {bb_d_D2, "Bethe-Bloch (Bohr)", kColBB, 1, true}, + {catima_d_D2, "CATIMA (default)", kColCATIMA, 2, true}, + {catimaPure_d_D2, "CATIMA (bare, no corr)", kColCATIMApure, 7, true}, + }; + if (table_d_D2) + ents.push_back({table_d_D2, "SRIM Table", kColTable, 3, false}); + DrawScenario("Deuteron in D_{2} (600 Torr)", ents, 0.2, 20.0, 10.0, c2); + } + + auto *c3 = new TCanvas("c3", "Alpha in D2", 1200, 900); + { + std::vector ents = { + {bb_a_D2, "Bethe-Bloch (Bohr)", kColBB, 1, true}, + {catima_a_D2, "CATIMA (default)", kColCATIMA, 2, true}, + {catimaPure_a_D2, "CATIMA (bare, no corr)", kColCATIMApure, 7, true}, + }; + if (table_a_D2) + ents.push_back({table_a_D2, "SRIM Table", kColTable, 3, false}); + DrawScenario("Alpha (^{4}He) in D_{2} (600 Torr)", ents, 0.5, 30.0, 10.0, c3); + } + + auto *c4 = new TCanvas("c4", "Proton in He", 1200, 900); + { + std::vector ents = { + {bb_p_He, "Bethe-Bloch (Bohr)", kColBB, 1, true}, + {catima_p_He, "CATIMA (default)", kColCATIMA, 2, true}, + {catimaPure_p_He, "CATIMA (bare, no corr)", kColCATIMApure, 7, true}, + }; + if (table_p_He) + ents.push_back({table_p_He, "SRIM Table", kColTable, 3, false}); + DrawScenario("Proton in ^{4}He (700 Torr)", ents, 0.1, 10.0, 5.0, c4); + } + + c1->Print("ELossComparison.pdf("); + c2->Print("ELossComparison.pdf"); + c3->Print("ELossComparison.pdf"); + c4->Print("ELossComparison.pdf)"); + + ::Info("ELossComparison", "Saved to ELossComparison.pdf"); +} diff --git a/macro/tests/ELossComparison/runTests.sh b/macro/tests/ELossComparison/runTests.sh new file mode 100755 index 000000000..04bf95d22 --- /dev/null +++ b/macro/tests/ELossComparison/runTests.sh @@ -0,0 +1,18 @@ +#!/bin/bash +# Run the ELossComparison integration macro. +# Must be called after sourcing build/config.sh (sets VMCWORKDIR, LD_LIBRARY_PATH, etc.) + +set -e + +SCRIPT_DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" +cd "${SCRIPT_DIR}" + +if [ -z "${VMCWORKDIR}" ]; then + echo "ERROR: VMCWORKDIR is not set. Source build/config.sh first." >&2 + exit 1 +fi + +echo "Running ELossComparison..." +root -l -b -q ELossComparison.C + +echo "Done. Output: ${SCRIPT_DIR}/ELossComparison.pdf" diff --git a/scripts/formatAll.sh b/scripts/formatAll.sh index 873e68e9d..9a9594886 100755 --- a/scripts/formatAll.sh +++ b/scripts/formatAll.sh @@ -2,16 +2,25 @@ #Break if something errors set -e -echo "Indenting all source (.cxx) and header files (.h)" +echo "Indenting all source (.cxx, .C) and header files (.h)" -#Make sure $VMCWORKDIR is set. +#Make sure $VMCWORKDIR is set. if [ -z "$VMCWORKDIR" ] then echo "VMCWORKDIR is unset. Please run config.sh on build directory first. Aborting..." return fi -find -L $VMCWORKDIR -type f | grep -w h$ | xargs -I fname clang-format -i --verbose fname -find -L $VMCWORKDIR -type f | grep -w cxx$ | xargs -I fname clang-format -i --verbose fname +# Directories excluded from formatting (matches .github/workflows/main.yml) +EXCLUDE_DIRS="macro scripts compiled resources include geometry build .venv lib bin" + +# Build a prune expression for find +PRUNE_EXPR="" +for dir in $EXCLUDE_DIRS; do + PRUNE_EXPR="$PRUNE_EXPR -path $VMCWORKDIR/$dir -prune -o" +done + +find -L $VMCWORKDIR $PRUNE_EXPR -type f \( -name "*.h" -o -name "*.cxx" -o -name "*.C" \) -print \ + | xargs -I fname clang-format-17 -i --verbose fname echo "Success! Indented all source and header files"