lars
62945c95c4
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.
87 lines
3.3 KiB
C++
87 lines
3.3 KiB
C++
// 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;
|
|
}
|