// 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;
}