forked from celeritas-project/celeritas
-
Notifications
You must be signed in to change notification settings - Fork 0
/
GeantMicroXsCalculator.cc
108 lines (96 loc) · 3.56 KB
/
GeantMicroXsCalculator.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
//----------------------------------*-C++-*----------------------------------//
// Copyright 2023-2024 UT-Battelle, LLC, and other Celeritas developers.
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: (Apache-2.0 OR MIT)
//---------------------------------------------------------------------------//
//! \file celeritas/ext/detail/GeantMicroXsCalculator.cc
//---------------------------------------------------------------------------//
#include "GeantMicroXsCalculator.hh"
#include <G4Material.hh>
#include <G4VEmModel.hh>
#include "corecel/cont/Range.hh"
#include "corecel/math/Algorithms.hh"
#include "celeritas/io/ImportUnits.hh"
namespace celeritas
{
namespace detail
{
//---------------------------------------------------------------------------//
/*!
* Construct with Geant4 classes.
*
* The secondary production cut should be nonzero if a secondary is produced.
*/
GeantMicroXsCalculator::GeantMicroXsCalculator(
G4VEmModel const& model,
G4ParticleDefinition const& particle,
G4Material const& material,
double secondary_production_cut)
: model_{const_cast<G4VEmModel&>(model)}
, particle_{particle}
, material_{material}
, secondary_cut_{secondary_production_cut}
{
CELER_EXPECT(secondary_cut_ >= 0);
}
//---------------------------------------------------------------------------//
/*!
* Calculate cross sections for all elements in the material.
*/
void GeantMicroXsCalculator::operator()(VecDouble const& energy_grid,
VecVecDouble* result_xs) const
{
CELER_EXPECT(result_xs);
auto const& elements = *material_.GetElementVector();
// Resize microscopic cross sections for all elements
result_xs->resize(elements.size());
for (auto& mxs_vec : *result_xs)
{
mxs_vec.resize(energy_grid.size());
}
auto calc_element_xs
= [this, &elements](std::size_t elcomp_idx, double energy) {
// Calculate microscopic cross section
double xs = model_.ComputeCrossSectionPerAtom(
&particle_,
elements[elcomp_idx],
energy,
secondary_cut_,
/* max_energy = */ std::numeric_limits<double>::max());
return clamp_to_nonneg(xs);
};
double const xs_scaling = native_value_from_clhep(ImportUnits::len_sq);
// Outer loop over energy to reduce material setup calls
for (auto energy_idx : range(energy_grid.size()))
{
double const energy = energy_grid[energy_idx];
CELER_ASSERT(energy > 0);
model_.SetupForMaterial(&particle_, &material_, energy);
// Inner loop over elements
for (auto elcomp_idx : range(elements.size()))
{
(*result_xs)[elcomp_idx][energy_idx]
= calc_element_xs(elcomp_idx, energy) * xs_scaling;
}
}
for (std::vector<double>& xs : *result_xs)
{
// Avoid cross-section vectors starting or ending with zero values.
// Geant4 simply uses the next/previous bin value when the vector's
// front/back are zero. This probably isn't correct but it replicates
// Geant4's behavior.
if (xs[0] == 0)
{
xs[0] = xs[1];
}
auto const last_idx = xs.size() - 1;
if (xs[last_idx] == 0)
{
// Cross-section ends with zero, use previous bin value
xs[last_idx] = xs[last_idx - 1];
}
}
}
//---------------------------------------------------------------------------//
} // namespace detail
} // namespace celeritas