From 62945c95c453a2e24525d5370bdd0a2a62531a34 Mon Sep 17 00:00:00 2001
From: Lars Bogner <upyxo@student.kit.edu>
Date: Thu, 25 Jun 2026 16:53:23 +0200
Subject: [PATCH] Add export_xsec executable for PbWO4 gamma attenuation
 coefficients

Builds a minimal PbWO4 geometry, forces EM physics tables via a
zero-event run, then sweeps photon energy 1 eV-10 MeV with
G4EmCalculator to dump per-process mass attenuation coefficients to CSV.
---
 CMakeLists.txt |  3 ++
 export_xsec.cc | 86 ++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 89 insertions(+)
 create mode 100644 export_xsec.cc

diff --git a/CMakeLists.txt b/CMakeLists.txt
index b6ce8ad..bbe7a56 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -63,6 +63,9 @@ target_link_libraries(run_pbwo4 ${Geant4_LIBRARIES} ${Python3_LIBRARIES})
 add_executable(run_sampling run_sampling.cc ${sources} ${headers})
 target_link_libraries(run_sampling ${Geant4_LIBRARIES} ${Python3_LIBRARIES})
 
+add_executable(export_xsec export_xsec.cc ${sources} ${headers})
+target_link_libraries(export_xsec ${Geant4_LIBRARIES} ${Python3_LIBRARIES})
+
 #----------------------------------------------------------------------------
 # Install the executable to 'bin' directory under CMAKE_INSTALL_PREFIX
 #
diff --git a/export_xsec.cc b/export_xsec.cc
new file mode 100644
index 0000000..3b45fbf
--- /dev/null
+++ b/export_xsec.cc
@@ -0,0 +1,86 @@
+// Export material-specific gamma cross sections for PbWO4 using G4EmCalculator.
+//
+// Builds a minimal geometry, initialises FTFP_BERT, forces the EM physics
+// tables to be built (via a fakeRun "/run/beamOn 0"), then sweeps photon
+// energy from 1 eV to 1e5 eV and writes a CSV of the mass attenuation
+// coefficients (cm^2/g) for each gamma process plus the total.
+
+#include "GeometryDescriptor.hh"
+#include "G4System.hh"
+
+#include "G4EmCalculator.hh"
+#include "G4Material.hh"
+#include "G4NistManager.hh"
+#include "G4Gamma.hh"
+#include "G4SystemOfUnits.hh"
+
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+int main(int argc, char** argv) {
+    const std::string materialName = "G4_PbWO4";
+    const std::string outfile = (argc > 1) ? argv[1] : "pbwo4_xsec.csv";
+
+    // Log-spaced energy grid: 1 eV -> 1e5 eV, points per decade (small steps).
+    const double eMin = 1.0 * eV;
+    const double eMax = 1.0e7 * eV;  // 10 MeV
+    const int pointsPerDecade = 200;
+    const int nDecades = 7;  // 1e0 .. 1e7 eV
+    const int nPoints = pointsPerDecade * nDecades + 1;
+
+    // Gamma processes registered by the EM standard physics in FTFP_BERT.
+    const std::vector<std::string> processes = {"phot", "compt", "conv", "Rayl"};
+
+    // Minimal valid geometry: one active PbWO4 layer so init() has something to build.
+    GeometryDescriptor gd;
+    gd.addLayer(2.0, materialName, true, 1, 1);
+
+    G4System g4;
+    g4.init(gd, -1);
+    // fakeRun: builds EM physics tables without invoking RunAction (no ROOT file).
+    g4.applyUICommand("/run/beamOn 0");
+
+    const G4Material* mat = G4NistManager::Instance()->FindOrBuildMaterial(materialName);
+    if (mat == nullptr) {
+        std::cerr << "Could not build material " << materialName << std::endl;
+        return 1;
+    }
+    const double density_g_cm3 = mat->GetDensity() / (g / cm3);
+
+    G4EmCalculator calc;
+    const G4ParticleDefinition* gamma = G4Gamma::Gamma();
+
+    std::ofstream out(outfile);
+    out << "# material=" << materialName
+        << " density_g_per_cm3=" << density_g_cm3 << "\n";
+    out << "# mass attenuation coefficients mu/rho in cm^2/g\n";
+    out << "energy_eV,phot,compt,conv,rayl,total\n";
+    out.setf(std::ios::scientific);
+    out.precision(6);
+
+    const double logStep = (std::log10(eMax) - std::log10(eMin)) / (nPoints - 1);
+    for (int i = 0; i < nPoints; ++i) {
+        const double energy = std::pow(10.0, std::log10(eMin) + i * logStep);
+
+        double total_lin = 0.0;  // linear attenuation, G4 units (1/mm)
+        std::vector<double> mu;  // per-process mass attenuation (cm^2/g)
+        mu.reserve(processes.size());
+        for (const auto& proc : processes) {
+            // Macroscopic cross section Sigma [1/mm] = material-specific cross section.
+            const double sigmaVol = calc.ComputeCrossSectionPerVolume(energy, gamma, proc, mat);
+            total_lin += sigmaVol;
+            // Convert to mass attenuation coefficient: (Sigma in 1/cm) / density.
+            mu.push_back((sigmaVol * cm) / density_g_cm3);
+        }
+        const double total_mu = (total_lin * cm) / density_g_cm3;
+
+        out << energy / eV << "," << mu[0] << "," << mu[1] << ","
+            << mu[2] << "," << mu[3] << "," << total_mu << "\n";
+    }
+
+    std::cout << "Wrote " << nPoints << " energy points for " << materialName
+              << " to " << outfile << std::endl;
+    return 0;
+}