• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

openmc-dev / openmc / 23016086442

12 Mar 2026 05:50PM UTC coverage: 81.564% (+0.001%) from 81.563%
23016086442

Pull #3874

github

web-flow
Merge ff9d6612d into 387b41ab6
Pull Request #3874: Fix a bug in compton scattering shell selection

17553 of 25264 branches covered (69.48%)

Branch coverage included in aggregate %.

28 of 29 new or added lines in 1 file covered. (96.55%)

1 existing line in 1 file now uncovered.

57960 of 67317 relevant lines covered (86.1%)

44967711.04 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

89.26
/src/photon.cpp
1
#include "openmc/photon.h"
2

3
#include "openmc/array.h"
4
#include "openmc/bremsstrahlung.h"
5
#include "openmc/constants.h"
6
#include "openmc/distribution_multi.h"
7
#include "openmc/hdf5_interface.h"
8
#include "openmc/message_passing.h"
9
#include "openmc/nuclide.h"
10
#include "openmc/particle.h"
11
#include "openmc/random_dist.h"
12
#include "openmc/random_lcg.h"
13
#include "openmc/search.h"
14
#include "openmc/settings.h"
15

16
#include "openmc/tensor.h"
17

18
#include <cmath>
19
#include <fmt/core.h>
20
#include <tuple> // for tie
21

22
namespace openmc {
23

24
constexpr int PhotonInteraction::MAX_STACK_SIZE;
25

26
//==============================================================================
27
// Global variables
28
//==============================================================================
29

30
namespace data {
31

32
tensor::Tensor<double> compton_profile_pz;
33

34
std::unordered_map<std::string, int> element_map;
35
vector<unique_ptr<PhotonInteraction>> elements;
36

37
} // namespace data
38

39
//==============================================================================
40
// PhotonInteraction implementation
41
//==============================================================================
42

43
PhotonInteraction::PhotonInteraction(hid_t group)
693✔
44
{
45
  // Set index of element in global vector
46
  index_ = data::elements.size();
693✔
47

48
  // Get name of nuclide from group, removing leading '/'
49
  name_ = object_name(group).substr(1);
693✔
50
  data::element_map[name_] = index_;
693✔
51

52
  // Get atomic number
53
  read_attribute(group, "Z", Z_);
693✔
54

55
  // Determine number of energies and read energy grid
56
  read_dataset(group, "energy", energy_);
693✔
57

58
  // Read coherent scattering
59
  hid_t rgroup = open_group(group, "coherent");
693✔
60
  read_dataset(rgroup, "xs", coherent_);
693✔
61

62
  hid_t dset = open_dataset(rgroup, "integrated_scattering_factor");
693✔
63
  coherent_int_form_factor_ = Tabulated1D {dset};
693✔
64
  close_dataset(dset);
693✔
65

66
  if (object_exists(group, "anomalous_real")) {
693!
67
    dset = open_dataset(rgroup, "anomalous_real");
×
68
    coherent_anomalous_real_ = Tabulated1D {dset};
×
69
    close_dataset(dset);
×
70
  }
71

72
  if (object_exists(group, "anomalous_imag")) {
693!
73
    dset = open_dataset(rgroup, "anomalous_imag");
×
74
    coherent_anomalous_imag_ = Tabulated1D {dset};
×
75
    close_dataset(dset);
×
76
  }
77
  close_group(rgroup);
693✔
78

79
  // Read incoherent scattering
80
  rgroup = open_group(group, "incoherent");
693✔
81
  read_dataset(rgroup, "xs", incoherent_);
693✔
82
  dset = open_dataset(rgroup, "scattering_factor");
693✔
83
  incoherent_form_factor_ = Tabulated1D {dset};
693✔
84
  close_dataset(dset);
693✔
85
  close_group(rgroup);
693✔
86

87
  // Read pair production
88
  if (object_exists(group, "pair_production_electron")) {
693!
89
    rgroup = open_group(group, "pair_production_electron");
693✔
90
    read_dataset(rgroup, "xs", pair_production_electron_);
693✔
91
    close_group(rgroup);
693✔
92
  } else {
93
    pair_production_electron_ = tensor::zeros_like(energy_);
×
94
  }
95

96
  // Read pair production
97
  if (object_exists(group, "pair_production_nuclear")) {
693!
98
    rgroup = open_group(group, "pair_production_nuclear");
693✔
99
    read_dataset(rgroup, "xs", pair_production_nuclear_);
693✔
100
    close_group(rgroup);
693✔
101
  } else {
102
    pair_production_nuclear_ = tensor::zeros_like(energy_);
×
103
  }
104

105
  // Read photoelectric
106
  rgroup = open_group(group, "photoelectric");
693✔
107
  read_dataset(rgroup, "xs", photoelectric_total_);
693✔
108
  close_group(rgroup);
693✔
109

110
  // Read heating
111
  if (object_exists(group, "heating")) {
693!
112
    rgroup = open_group(group, "heating");
693✔
113
    read_dataset(rgroup, "xs", heating_);
693✔
114
    close_group(rgroup);
693✔
115
  } else {
116
    heating_ = tensor::zeros_like(energy_);
×
117
  }
118

119
  // Read subshell photoionization cross section and atomic relaxation data
120
  rgroup = open_group(group, "subshells");
693✔
121
  vector<std::string> designators;
693✔
122
  read_attribute(rgroup, "designators", designators);
693✔
123
  auto n_shell = designators.size();
693!
124
  if (n_shell == 0) {
693!
125
    throw std::runtime_error {
×
126
      "Photoatomic data for " + name_ + " does not have subshell data."};
×
127
  }
128

129
  shells_.resize(n_shell);
693✔
130
  cross_sections_ = tensor::zeros<double>({energy_.size(), n_shell});
693✔
131

132
  // Create mapping from designator to index
133
  std::unordered_map<int, int> shell_map;
693✔
134
  for (int i = 0; i < n_shell; ++i) {
6,232✔
135
    const auto& designator {designators[i]};
5,539✔
136

137
    int j = 1;
5,539✔
138
    for (const auto& subshell : SUBSHELLS) {
41,609!
139
      if (designator == subshell) {
41,609✔
140
        shell_map[j] = i;
5,539✔
141
        shells_[i].index_subshell = j;
5,539✔
142
        break;
5,539✔
143
      }
144
      ++j;
36,070✔
145
    }
146
  }
147
  shell_map[0] = -1;
693✔
148

149
  for (int i = 0; i < n_shell; ++i) {
6,232✔
150
    const auto& designator {designators[i]};
5,539✔
151
    auto& shell {shells_[i]};
5,539✔
152

153
    // TODO: Move to ElectronSubshell constructor
154

155
    hid_t tgroup = open_group(rgroup, designator.c_str());
5,539✔
156

157
    // Read binding energy if atomic relaxation data is present
158
    if (attribute_exists(tgroup, "binding_energy")) {
5,539!
159
      has_atomic_relaxation_ = true;
5,539✔
160
      read_attribute(tgroup, "binding_energy", shell.binding_energy);
5,539✔
161
    }
162

163
    // Read subshell cross section
164
    tensor::Tensor<double> xs;
5,539✔
165
    dset = open_dataset(tgroup, "xs");
5,539✔
166
    read_attribute(dset, "threshold_idx", shell.threshold);
5,539✔
167
    close_dataset(dset);
5,539✔
168
    read_dataset(tgroup, "xs", xs);
5,539✔
169

170
    auto cross_section =
5,539✔
171
      cross_sections_.slice(tensor::range(static_cast<size_t>(shell.threshold),
5,539✔
172
                              cross_sections_.shape(0)),
173
        i);
11,078!
174
    cross_section = tensor::where(xs > 0, tensor::log(xs), 0);
16,617✔
175

176
    if (settings::atomic_relaxation && object_exists(tgroup, "transitions")) {
5,539✔
177
      // Determine dimensions of transitions
178
      dset = open_dataset(tgroup, "transitions");
3,139✔
179
      auto dims = object_shape(dset);
3,139✔
180
      close_dataset(dset);
3,139✔
181

182
      int n_transition = dims[0];
3,139!
183
      if (n_transition > 0) {
3,139!
184
        tensor::Tensor<double> matrix;
3,139✔
185
        read_dataset(tgroup, "transitions", matrix);
3,139✔
186

187
        // Transition probability normalization
188
        double norm =
3,139✔
189
          tensor::Tensor<double>(matrix.slice(tensor::all, 3)).sum();
6,278✔
190

191
        shell.transitions.resize(n_transition);
3,139✔
192
        for (int j = 0; j < n_transition; ++j) {
167,648✔
193
          auto& transition = shell.transitions[j];
164,509✔
194
          transition.primary_subshell = shell_map.at(matrix(j, 0));
164,509✔
195
          transition.secondary_subshell = shell_map.at(matrix(j, 1));
164,509✔
196
          transition.energy = matrix(j, 2);
164,509✔
197
          transition.probability = matrix(j, 3) / norm;
164,509✔
198
        }
199
      }
3,139✔
200
    }
3,139✔
201
    close_group(tgroup);
5,539✔
202
  }
11,078✔
203
  close_group(rgroup);
693✔
204

205
  // Check the maximum size of the atomic relaxation stack
206
  auto max_size = this->calc_max_stack_size();
693✔
207
  if (max_size > MAX_STACK_SIZE && mpi::master) {
693!
208
    warning(fmt::format("The subshell vacancy stack in atomic relaxation can "
×
209
                        "grow up to {}, but the stack size limit is set to {}.",
210
      max_size, MAX_STACK_SIZE));
211
  }
212

213
  // Determine number of electron shells
214
  rgroup = open_group(group, "compton_profiles");
693✔
215

216
  // Read electron shell PDF and binding energies
217
  read_dataset(rgroup, "num_electrons", electron_pdf_);
693✔
218
  electron_pdf_ /= electron_pdf_.sum();
693✔
219
  read_dataset(rgroup, "binding_energy", binding_energy_);
693✔
220

221
  // Read Compton profiles
222
  read_dataset(rgroup, "J", profile_pdf_);
693✔
223

224
  // Get Compton profile momentum grid
225
  if (data::compton_profile_pz.size() == 0) {
693✔
226
    read_dataset(rgroup, "pz", data::compton_profile_pz);
312✔
227
  }
228
  close_group(rgroup);
693✔
229

230
  // Map Compton subshell data to atomic relaxation data by finding the
231
  // subshell with the equivalent binding energy
232
  if (settings::atomic_relaxation && has_atomic_relaxation_) {
693!
233
    auto is_close = [](double a, double b) {
9,652✔
234
      return std::abs(a - b) / a < FP_REL_PRECISION;
8,974✔
235
    };
236
    subshell_map_ = tensor::Tensor<int>(binding_energy_.shape(), -1);
1,356✔
237
    for (int i = 0; i < binding_energy_.size(); ++i) {
4,987✔
238
      double E_b = binding_energy_[i];
4,309!
239
      if (i < n_shell && is_close(E_b, shells_[i].binding_energy)) {
4,309!
240
        subshell_map_[i] = i;
3,570✔
241
      } else {
242
        for (int j = 0; j < n_shell; ++j) {
4,665!
243
          if (is_close(E_b, shells_[j].binding_energy)) {
4,665✔
244
            subshell_map_[i] = j;
739✔
245
            break;
739✔
246
          }
247
        }
248
      }
249
    }
250
  }
251

252
  // Create Compton profile CDF
253
  auto n_profile = data::compton_profile_pz.size();
693!
254
  auto n_shell_compton = profile_pdf_.shape(0);
693!
255
  profile_cdf_ = tensor::Tensor<double>({n_shell_compton, n_profile});
693✔
256
  for (int i = 0; i < n_shell_compton; ++i) {
5,347✔
257
    double c = 0.0;
4,654✔
258
    profile_cdf_(i, 0) = 0.0;
4,654✔
259
    for (int j = 0; j < n_profile - 1; ++j) {
144,274✔
260
      c += 0.5 *
139,620✔
261
           (data::compton_profile_pz(j + 1) - data::compton_profile_pz(j)) *
139,620✔
262
           (profile_pdf_(i, j) + profile_pdf_(i, j + 1));
139,620✔
263
      profile_cdf_(i, j + 1) = c;
139,620✔
264
    }
265
  }
266

267
  // Calculate total pair production
268
  pair_production_total_ = pair_production_nuclear_ + pair_production_electron_;
693✔
269

270
  if (settings::electron_treatment == ElectronTreatment::TTB) {
693✔
271
    // Read bremsstrahlung scaled DCS
272
    rgroup = open_group(group, "bremsstrahlung");
557✔
273
    read_dataset(rgroup, "dcs", dcs_);
557✔
274
    auto n_e = dcs_.shape(0);
557!
275
    auto n_k = dcs_.shape(1);
557!
276

277
    // Get energy grids used for bremsstrahlung DCS and for stopping powers
278
    tensor::Tensor<double> electron_energy;
557✔
279
    read_dataset(rgroup, "electron_energy", electron_energy);
557✔
280
    if (data::ttb_k_grid.size() == 0) {
557✔
281
      read_dataset(rgroup, "photon_energy", data::ttb_k_grid);
264✔
282
    }
283

284
    // Get data used for density effect correction
285
    read_dataset(rgroup, "num_electrons", n_electrons_);
557✔
286
    read_dataset(rgroup, "ionization_energy", ionization_energy_);
557✔
287
    read_attribute(rgroup, "I", I_);
557✔
288
    close_group(rgroup);
557✔
289

290
    // Truncate the bremsstrahlung data at the cutoff energy
291
    int photon = ParticleType::photon().transport_index();
557✔
292
    const auto& E {electron_energy};
557✔
293
    double cutoff = settings::energy_cutoff[photon];
557✔
294
    if (cutoff > E(0)) {
557✔
295
      size_t i_grid = lower_bound_index(
11✔
296
        E.cbegin(), E.cend(), settings::energy_cutoff[photon]);
11✔
297

298
      // calculate interpolation factor
299
      double f = (std::log(cutoff) - std::log(E(i_grid))) /
11✔
300
                 (std::log(E(i_grid + 1)) - std::log(E(i_grid)));
11✔
301

302
      // Interpolate bremsstrahlung DCS at the cutoff energy and truncate
303
      tensor::Tensor<double> dcs({n_e - i_grid, n_k});
11✔
304
      for (int i = 0; i < n_k; ++i) {
341✔
305
        double y = std::exp(
990✔
306
          std::log(dcs_(i_grid, i)) +
330✔
307
          f * (std::log(dcs_(i_grid + 1, i)) - std::log(dcs_(i_grid, i))));
330✔
308
        tensor::View<double> col_i = dcs.slice(tensor::all, i);
330✔
309
        col_i(0) = y;
330✔
310
        for (int j = i_grid + 1; j < n_e; ++j) {
61,290✔
311
          col_i(j - i_grid) = dcs_(j, i);
60,960✔
312
        }
313
      }
330✔
314
      dcs_ = dcs;
11✔
315

316
      tensor::Tensor<double> frst({static_cast<size_t>(1)});
11✔
317
      frst(0) = cutoff;
11✔
318
      tensor::Tensor<double> rest(electron_energy.slice(
11✔
319
        tensor::range(i_grid + 1, electron_energy.size())));
11✔
320
      electron_energy = tensor::concatenate(frst, rest);
22✔
321
    }
33✔
322

323
    // Set incident particle energy grid
324
    if (data::ttb_e_grid.size() == 0) {
557✔
325
      data::ttb_e_grid = electron_energy;
264✔
326
    }
327

328
    // Calculate the radiative stopping power
329
    stopping_power_radiative_ =
557✔
330
      tensor::Tensor<double>({data::ttb_e_grid.size()});
557✔
331
    for (int i = 0; i < data::ttb_e_grid.size(); ++i) {
111,800✔
332
      // Integrate over reduced photon energy
333
      double c = 0.0;
334
      for (int j = 0; j < data::ttb_k_grid.size() - 1; ++j) {
3,337,290✔
335
        c += 0.5 * (dcs_(i, j + 1) + dcs_(i, j)) *
3,226,047✔
336
             (data::ttb_k_grid(j + 1) - data::ttb_k_grid(j));
3,226,047✔
337
      }
338
      double e = data::ttb_e_grid(i);
111,243✔
339

340
      // Square of the ratio of the speed of light to the velocity of the
341
      // charged particle
342
      double beta_sq = e * (e + 2.0 * MASS_ELECTRON_EV) /
111,243✔
343
                       ((e + MASS_ELECTRON_EV) * (e + MASS_ELECTRON_EV));
111,243✔
344

345
      stopping_power_radiative_(i) = Z_ * Z_ / beta_sq * e * c;
111,243✔
346
    }
347
  }
557✔
348

349
  // Take logarithm of energies and cross sections since they are log-log
350
  // interpolated. Note that cross section libraries converted from ACE files
351
  // represent zero as exp(-500) to avoid log-log interpolation errors. For
352
  // values below exp(-499) we store the log as -900, for which exp(-900)
353
  // evaluates to zero.
354
  double limit = std::exp(-499.0);
693✔
355
  energy_ = tensor::log(energy_);
693✔
356
  coherent_ = tensor::where(coherent_ > limit, tensor::log(coherent_), -900.0);
2,079✔
357
  incoherent_ =
693✔
358
    tensor::where(incoherent_ > limit, tensor::log(incoherent_), -900.0);
2,079✔
359
  photoelectric_total_ = tensor::where(
1,386✔
360
    photoelectric_total_ > limit, tensor::log(photoelectric_total_), -900.0);
2,079✔
361
  pair_production_total_ = tensor::where(pair_production_total_ > limit,
2,079✔
362
    tensor::log(pair_production_total_), -900.0);
1,386✔
363
  heating_ = tensor::where(heating_ > limit, tensor::log(heating_), -900.0);
2,772✔
364
}
693✔
365

366
PhotonInteraction::~PhotonInteraction()
693✔
367
{
368
  data::element_map.erase(name_);
693✔
369
}
13,167✔
370

371
int PhotonInteraction::calc_max_stack_size() const
693✔
372
{
373
  // Table to store solutions to sub-problems
374
  std::unordered_map<int, int> visited;
693✔
375

376
  // Find the maximum possible size of the stack used to store holes created
377
  // during atomic relaxation, checking over every subshell the initial hole
378
  // could be in
379
  int max_size = 0;
693✔
380
  for (int i_shell = 0; i_shell < shells_.size(); ++i_shell) {
6,232✔
381
    max_size = std::max(max_size, this->calc_helper(visited, i_shell));
6,232✔
382
  }
383
  return max_size;
693✔
384
}
693✔
385

386
int PhotonInteraction::calc_helper(
319,029✔
387
  std::unordered_map<int, int>& visited, int i_shell) const
388
{
389
  // No transitions for this subshell, so this is the only shell in the stack
390
  const auto& shell {shells_[i_shell]};
319,029✔
391
  if (shell.transitions.empty()) {
319,029✔
392
    return 1;
393
  }
394

395
  // Check the table to see if the maximum stack size has already been
396
  // calculated for this shell
397
  auto it = visited.find(i_shell);
187,636✔
398
  if (it != visited.end()) {
187,636✔
399
    return it->second;
184,497✔
400
  }
401

402
  int max_size = 0;
3,139✔
403
  for (const auto& transition : shell.transitions) {
167,648✔
404
    // If this is a non-radiative transition two vacancies are created and
405
    // the stack grows by one; if this is a radiative transition only one
406
    // vacancy is created and the stack size stays the same
407
    int size = 0;
164,509✔
408
    if (transition.secondary_subshell != -1) {
164,509✔
409
      size = this->calc_helper(visited, transition.secondary_subshell) + 1;
148,981✔
410
    }
411
    size =
329,018✔
412
      std::max(size, this->calc_helper(visited, transition.primary_subshell));
164,509✔
413
    max_size = std::max(max_size, size);
169,479✔
414
  }
415
  visited[i_shell] = max_size;
3,139✔
416
  return max_size;
3,139✔
417
}
418

419
void PhotonInteraction::compton_scatter(double alpha, bool doppler,
19,583,056✔
420
  double* alpha_out, double* mu, int* i_shell, uint64_t* seed) const
421
{
422
  double form_factor_xmax = 0.0;
19,583,056✔
423
  while (true) {
20,057,844✔
424
    // Sample Klein-Nishina distribution for trial energy and angle
425
    std::tie(*alpha_out, *mu) = klein_nishina(alpha, seed);
20,057,844✔
426

427
    // Note that the parameter used here does not correspond exactly to the
428
    // momentum transfer q in ENDF-102 Eq. (27.2). Rather, this is the
429
    // parameter as defined by Hubbell, where the actual data comes from
430
    double x =
20,057,844✔
431
      MASS_ELECTRON_EV / PLANCK_C * alpha * std::sqrt(0.5 * (1.0 - *mu));
20,057,844✔
432

433
    // Calculate S(x, Z) and S(x_max, Z)
434
    double form_factor_x = incoherent_form_factor_(x);
20,057,844✔
435
    if (form_factor_xmax == 0.0) {
20,057,844✔
436
      form_factor_xmax =
19,583,056✔
437
        incoherent_form_factor_(MASS_ELECTRON_EV / PLANCK_C * alpha);
19,583,056✔
438
    }
439

440
    // Perform rejection on form factor
441
    if (prn(seed) < form_factor_x / form_factor_xmax) {
20,057,844✔
442
      if (doppler) {
19,583,056!
443
        double E_out;
19,583,056✔
444
        this->compton_doppler(alpha, *mu, &E_out, i_shell, seed);
19,583,056✔
445
        *alpha_out = E_out / MASS_ELECTRON_EV;
19,583,056✔
446
      } else {
NEW
447
        *i_shell = -1;
×
448
      }
449
      break;
19,583,056✔
450
    }
451
  }
452
}
19,583,056✔
453

454
void PhotonInteraction::compton_doppler(
19,583,056✔
455
  double alpha, double mu, double* E_out, int* i_shell, uint64_t* seed) const
456
{
457
  auto n = data::compton_profile_pz.size();
19,583,056✔
458

459
  int shell; // index for shell
19,584,266✔
460
  while (true) {
19,584,266✔
461
    // Sample electron shell
462
    double rn = prn(seed);
19,584,266✔
463
    double c = 0.0;
19,584,266✔
464
    for (shell = 0; shell < electron_pdf_.size(); ++shell) {
52,195,034!
465
      c += electron_pdf_(shell);
52,195,034✔
466
      if (rn < c)
52,195,034✔
467
        break;
468
    }
469

470
    // Determine binding energy of shell
471
    double E_b = binding_energy_(shell);
19,584,266✔
472

473
    // Determine p_z,max
474
    double E = alpha * MASS_ELECTRON_EV;
19,584,266✔
475
    if (E < E_b) {
19,584,266✔
476
      continue;
1,210✔
477
    }
478

479
    double pz_max = -FINE_STRUCTURE * (E_b - (E - E_b) * alpha * (1.0 - mu)) /
39,166,112✔
480
                    std::sqrt(2.0 * E * (E - E_b) * (1.0 - mu) + E_b * E_b);
19,583,056✔
481
    if (pz_max < 0.0) {
19,583,056✔
482
      *E_out = alpha / (1 + alpha * (1 - mu)) * MASS_ELECTRON_EV;
182,754✔
483
      break;
182,754✔
484
    }
485

486
    // Determine profile cdf value corresponding to p_z,max
487
    double c_max;
19,400,302✔
488
    if (pz_max > data::compton_profile_pz(n - 1)) {
19,400,302✔
489
      c_max = profile_cdf_(shell, n - 1);
1,370,667✔
490
    } else {
491
      int i = lower_bound_index(data::compton_profile_pz.cbegin(),
18,029,635✔
492
        data::compton_profile_pz.cend(), pz_max);
18,029,635✔
493
      double pz_l = data::compton_profile_pz(i);
18,029,635!
494
      double pz_r = data::compton_profile_pz(i + 1);
18,029,635✔
495
      double p_l = profile_pdf_(shell, i);
18,029,635✔
496
      double p_r = profile_pdf_(shell, i + 1);
18,029,635✔
497
      double c_l = profile_cdf_(shell, i);
18,029,635✔
498
      if (pz_l == pz_r) {
18,029,635!
499
        c_max = c_l;
500
      } else if (p_l == p_r) {
18,029,635✔
501
        c_max = c_l + (pz_max - pz_l) * p_l;
26,964✔
502
      } else {
503
        double m = (p_l - p_r) / (pz_l - pz_r);
18,002,671✔
504
        c_max = c_l + (std::pow((m * (pz_max - pz_l) + p_l), 2) - p_l * p_l) /
18,002,671✔
505
                        (2.0 * m);
18,002,671✔
506
      }
507
    }
508

509
    // Sample value on bounded cdf
510
    c = prn(seed) * c_max;
19,400,302✔
511

512
    // Determine pz corresponding to sampled cdf value
513
    tensor::View<const double> cdf_shell = profile_cdf_.slice(shell);
19,400,302✔
514
    int i = lower_bound_index(cdf_shell.cbegin(), cdf_shell.cend(), c);
19,400,302✔
515
    double pz_l = data::compton_profile_pz(i);
19,400,302!
516
    double pz_r = data::compton_profile_pz(i + 1);
19,400,302✔
517
    double p_l = profile_pdf_(shell, i);
19,400,302✔
518
    double p_r = profile_pdf_(shell, i + 1);
19,400,302✔
519
    double c_l = profile_cdf_(shell, i);
19,400,302✔
520
    double pz;
19,400,302✔
521
    if (pz_l == pz_r) {
19,400,302!
522
      pz = pz_l;
523
    } else if (p_l == p_r) {
19,400,302✔
524
      pz = pz_l + (c - c_l) / p_l;
1,521,581✔
525
    } else {
526
      double m = (p_l - p_r) / (pz_l - pz_r);
17,878,721✔
527
      pz = pz_l + (std::sqrt(p_l * p_l + 2.0 * m * (c - c_l)) - p_l) / m;
17,878,721✔
528
    }
529

530
    // Determine outgoing photon energy corresponding to electron momentum
531
    // (solve Eq. 39 in LA-UR-04-0487 for E')
532
    double momentum_sq = std::pow((pz / FINE_STRUCTURE), 2);
19,400,302✔
533
    double f = 1.0 + alpha * (1.0 - mu);
19,400,302✔
534
    double a = momentum_sq - f * f;
19,400,302✔
535
    double b = 2.0 * E * (f - momentum_sq * mu);
19,400,302✔
536
    c = E * E * (momentum_sq - 1.0);
19,400,302✔
537

538
    double quad = b * b - 4.0 * a * c;
19,400,302✔
539
    if (quad < 0) {
19,400,302!
540
      *E_out = alpha / (1 + alpha * (1 - mu)) * MASS_ELECTRON_EV;
×
541
      break;
×
542
    }
543
    quad = std::sqrt(quad);
19,400,302✔
544
    double E_out1 = -(b + quad) / (2.0 * a);
19,400,302✔
545
    double E_out2 = -(b - quad) / (2.0 * a);
19,400,302✔
546

547
    // Determine solution to quadratic equation that is positive
548
    if (E_out1 > 0.0) {
19,400,302!
549
      if (E_out2 > 0.0) {
19,400,302!
550
        // If both are positive, pick one at random
551
        *E_out = prn(seed) < 0.5 ? E_out1 : E_out2;
29,108,040✔
552
      } else {
553
        *E_out = E_out1;
×
554
      }
555
    } else {
556
      if (E_out2 > 0.0) {
×
557
        *E_out = E_out2;
×
558
      } else {
559
        // No positive solution -- resample
560
        continue;
×
561
      }
562
    }
563
    if (*E_out < E - E_b)
19,400,302!
564
      break;
565
  }
1,210✔
566

567
  *i_shell = shell;
19,583,056✔
568
}
19,583,056✔
569

570
void PhotonInteraction::calculate_xs(Particle& p) const
60,957,353✔
571
{
572
  // Perform binary search on the element energy grid in order to determine
573
  // which points to interpolate between
574
  int n_grid = energy_.size();
60,957,353✔
575
  double log_E = std::log(p.E());
60,957,353✔
576
  int i_grid;
60,957,353✔
577
  if (log_E <= energy_[0]) {
60,957,353!
578
    i_grid = 0;
579
  } else if (log_E > energy_(n_grid - 1)) {
60,957,353!
580
    i_grid = n_grid - 2;
×
581
  } else {
582
    // We use upper_bound_index here because sometimes photons are created with
583
    // energies that exactly match a grid point
584
    i_grid = upper_bound_index(energy_.cbegin(), energy_.cend(), log_E);
60,957,353✔
585
  }
586

587
  // check for case where two energy points are the same
588
  if (energy_(i_grid) == energy_(i_grid + 1))
60,957,353!
589
    ++i_grid;
×
590

591
  // calculate interpolation factor
592
  double f =
60,957,353✔
593
    (log_E - energy_(i_grid)) / (energy_(i_grid + 1) - energy_(i_grid));
60,957,353✔
594

595
  auto& xs {p.photon_xs(index_)};
60,957,353✔
596
  xs.index_grid = i_grid;
60,957,353✔
597
  xs.interp_factor = f;
60,957,353✔
598

599
  // Calculate microscopic coherent cross section
600
  xs.coherent = std::exp(
121,914,706✔
601
    coherent_(i_grid) + f * (coherent_(i_grid + 1) - coherent_(i_grid)));
60,957,353✔
602

603
  // Calculate microscopic incoherent cross section
604
  xs.incoherent = std::exp(
121,914,706✔
605
    incoherent_(i_grid) + f * (incoherent_(i_grid + 1) - incoherent_(i_grid)));
60,957,353✔
606

607
  // Calculate microscopic photoelectric cross section
608
  xs.photoelectric = 0.0;
60,957,353✔
609
  tensor::View<const double> xs_lower = cross_sections_.slice(i_grid);
60,957,353✔
610
  tensor::View<const double> xs_upper = cross_sections_.slice(i_grid + 1);
60,957,353✔
611

612
  for (int i = 0; i < xs_upper.size(); ++i)
805,573,162✔
613
    if (xs_lower(i) != 0)
341,829,228✔
614
      xs.photoelectric +=
331,377,252✔
615
        std::exp(xs_lower(i) + f * (xs_upper(i) - xs_lower(i)));
331,377,252✔
616

617
  // Calculate microscopic pair production cross section
618
  xs.pair_production = std::exp(
121,914,706✔
619
    pair_production_total_(i_grid) +
60,957,353✔
620
    f * (pair_production_total_(i_grid + 1) - pair_production_total_(i_grid)));
60,957,353✔
621

622
  // Calculate microscopic total cross section
623
  xs.total =
60,957,353✔
624
    xs.coherent + xs.incoherent + xs.photoelectric + xs.pair_production;
60,957,353✔
625
  xs.last_E = p.E();
60,957,353✔
626
}
121,914,706✔
627

628
double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t* seed) const
1,555,934✔
629
{
630
  double mu;
1,780,037✔
631
  while (true) {
2,004,140✔
632
    // Determine maximum value of x^2
633
    double x2_max = std::pow(MASS_ELECTRON_EV / PLANCK_C * alpha, 2);
1,780,037✔
634

635
    // Determine F(x^2_max, Z)
636
    double F_max = coherent_int_form_factor_(x2_max);
1,780,037✔
637

638
    // Sample cumulative distribution
639
    double F = prn(seed) * F_max;
1,780,037✔
640

641
    // Determine x^2 corresponding to F
642
    const auto& x {coherent_int_form_factor_.x()};
1,780,037✔
643
    const auto& y {coherent_int_form_factor_.y()};
1,780,037✔
644
    int i = lower_bound_index(y.cbegin(), y.cend(), F);
1,780,037✔
645
    double r = (F - y[i]) / (y[i + 1] - y[i]);
1,780,037✔
646
    double x2 = x[i] + r * (x[i + 1] - x[i]);
1,780,037✔
647

648
    // Calculate mu
649
    mu = 1.0 - 2.0 * x2 / x2_max;
1,780,037✔
650

651
    if (prn(seed) < 0.5 * (1.0 + mu * mu))
1,780,037✔
652
      break;
653
  }
224,103✔
654
  return mu;
1,555,934✔
655
}
656

657
void PhotonInteraction::pair_production(double alpha, double* E_electron,
92,661✔
658
  double* E_positron, double* mu_electron, double* mu_positron,
659
  uint64_t* seed) const
660
{
661
  constexpr double r[] {122.81, 73.167, 69.228, 67.301, 64.696, 61.228, 57.524,
92,661✔
662
    54.033, 50.787, 47.851, 46.373, 45.401, 44.503, 43.815, 43.074, 42.321,
663
    41.586, 40.953, 40.524, 40.256, 39.756, 39.144, 38.462, 37.778, 37.174,
664
    36.663, 35.986, 35.317, 34.688, 34.197, 33.786, 33.422, 33.068, 32.740,
665
    32.438, 32.143, 31.884, 31.622, 31.438, 31.142, 30.950, 30.758, 30.561,
666
    30.285, 30.097, 29.832, 29.581, 29.411, 29.247, 29.085, 28.930, 28.721,
667
    28.580, 28.442, 28.312, 28.139, 27.973, 27.819, 27.675, 27.496, 27.285,
668
    27.093, 26.911, 26.705, 26.516, 26.304, 26.108, 25.929, 25.730, 25.577,
669
    25.403, 25.245, 25.100, 24.941, 24.790, 24.655, 24.506, 24.391, 24.262,
670
    24.145, 24.039, 23.922, 23.813, 23.712, 23.621, 23.523, 23.430, 23.331,
671
    23.238, 23.139, 23.048, 22.967, 22.833, 22.694, 22.624, 22.545, 22.446,
672
    22.358, 22.264};
673

674
  // The reduced screening radius r is the ratio of the screening radius to
675
  // the Compton wavelength of the electron, where the screening radius is
676
  // obtained under the assumption that the Coulomb field of the nucleus is
677
  // exponentially screened by atomic electrons. This allows us to use a
678
  // simplified atomic form factor and analytical approximations of the
679
  // screening functions in the pair production DCS instead of computing the
680
  // screening functions numerically. The reduced screening radii above for
681
  // Z = 1-99 come from F. Salvat, J. M. Fernández-Varea, and J. Sempau,
682
  // "PENELOPE-2011: A Code System for Monte Carlo Simulation of Electron and
683
  // Photon Transport," OECD-NEA, Issy-les-Moulineaux, France (2011).
684

685
  // Compute the high-energy Coulomb correction
686
  double a = Z_ / FINE_STRUCTURE;
92,661✔
687
  double c =
92,661✔
688
    a * a *
92,661✔
689
    (1.0 / (1.0 + a * a) + 0.202059 +
92,661✔
690
      a * a *
92,661✔
691
        (-0.03693 +
92,661✔
692
          a * a *
92,661✔
693
            (0.00835 +
92,661✔
694
              a * a *
92,661✔
695
                (-0.00201 +
92,661✔
696
                  a * a * (0.00049 + a * a * (-0.00012 + a * a * 0.00003))))));
92,661✔
697

698
  // The analytical approximation of the DCS underestimates the cross section
699
  // at low energies. The correction factor f compensates for this.
700
  double q = std::sqrt(2.0 / alpha);
92,661✔
701
  double f = q * (-0.1774 - 12.10 * a + 11.18 * a * a) +
92,661✔
702
             q * q * (8.523 + 73.26 * a - 44.41 * a * a) +
92,661✔
703
             q * q * q * (-13.52 - 121.1 * a + 96.41 * a * a) +
92,661✔
704
             q * q * q * q * (8.946 + 62.05 * a - 63.41 * a * a);
92,661✔
705

706
  // Calculate phi_1(1/2) and phi_2(1/2). The unnormalized PDF for the reduced
707
  // energy is given by p = 2*(1/2 - e)^2*phi_1(e) + phi_2(e), where phi_1 and
708
  // phi_2 are non-negative and maximum at e = 1/2.
709
  double b = 2.0 * r[Z_] / alpha;
92,661✔
710
  double t1 = 2.0 * std::log(1.0 + b * b);
92,661✔
711
  double t2 = b * std::atan(1.0 / b);
92,661✔
712
  double t3 = b * b * (4.0 - 4.0 * t2 - 3.0 * std::log(1.0 + 1.0 / (b * b)));
92,661✔
713
  double t4 = 4.0 * std::log(r[Z_]) - 4.0 * c + f;
92,661✔
714
  double phi1_max = 7.0 / 3.0 - t1 - 6.0 * t2 - t3 + t4;
92,661✔
715
  double phi2_max = 11.0 / 6.0 - t1 - 3.0 * t2 + 0.5 * t3 + t4;
92,661✔
716

717
  // To aid sampling, the unnormalized PDF can be expressed as
718
  // p = u_1*U_1(e)*pi_1(e) + u_2*U_2(e)*pi_2(e), where pi_1 and pi_2 are
719
  // normalized PDFs on the interval (e_min, e_max) from which values of e can
720
  // be sampled using the inverse transform method, and
721
  // U_1 = phi_1(e)/phi_1(1/2) and U_2 = phi_2(e)/phi_2(1/2) are valid
722
  // rejection functions. The reduced energy can now be sampled using a
723
  // combination of the composition and rejection methods.
724
  double u1 = 2.0 / 3.0 * std::pow(0.5 - 1.0 / alpha, 2) * phi1_max;
92,661✔
725
  double u2 = phi2_max;
92,661✔
726
  double e;
116,333✔
727
  while (true) {
116,333✔
728
    double rn = prn(seed);
116,333✔
729

730
    // Sample the index i in (1, 2) using the point probabilities
731
    // p(1) = u_1/(u_1 + u_2) and p(2) = u_2/(u_1 + u_2)
732
    int i;
116,333✔
733
    if (prn(seed) < u1 / (u1 + u2)) {
116,333✔
734
      i = 1;
10,428✔
735

736
      // Sample e from pi_1 using the inverse transform method
737
      e = rn >= 0.5
10,428✔
738
            ? 0.5 + (0.5 - 1.0 / alpha) * std::pow(2.0 * rn - 1.0, 1.0 / 3.0)
10,428✔
739
            : 0.5 - (0.5 - 1.0 / alpha) * std::pow(1.0 - 2.0 * rn, 1.0 / 3.0);
5,335✔
740
    } else {
741
      i = 2;
105,905✔
742

743
      // Sample e from pi_2 using the inverse transform method
744
      e = 1.0 / alpha + (0.5 - 1.0 / alpha) * 2.0 * rn;
105,905✔
745
    }
746

747
    // Calculate phi_i(e) and deliver e if rn <= U_i(e)
748
    b = r[Z_] / (2.0 * alpha * e * (1.0 - e));
116,333✔
749
    t1 = 2.0 * std::log(1.0 + b * b);
116,333✔
750
    t2 = b * std::atan(1.0 / b);
116,333✔
751
    t3 = b * b * (4.0 - 4.0 * t2 - 3.0 * std::log(1.0 + 1.0 / (b * b)));
116,333✔
752
    if (i == 1) {
116,333✔
753
      double phi1 = 7.0 / 3.0 - t1 - 6.0 * t2 - t3 + t4;
10,428✔
754
      if (prn(seed) <= phi1 / phi1_max)
10,428✔
755
        break;
756
    } else {
757
      double phi2 = 11.0 / 6.0 - t1 - 3.0 * t2 + 0.5 * t3 + t4;
105,905✔
758
      if (prn(seed) <= phi2 / phi2_max)
105,905✔
759
        break;
760
    }
761
  }
762

763
  // Compute the kinetic energy of the electron and the positron
764
  *E_electron = (alpha * e - 1.0) * MASS_ELECTRON_EV;
92,661✔
765
  *E_positron = (alpha * (1.0 - e) - 1.0) * MASS_ELECTRON_EV;
92,661✔
766

767
  // Sample the scattering angle of the electron. The cosine of the polar
768
  // angle of the direction relative to the incident photon is sampled from
769
  // p(mu) = C/(1 - beta*mu)^2 using the inverse transform method.
770
  double beta =
92,661✔
771
    std::sqrt(*E_electron * (*E_electron + 2.0 * MASS_ELECTRON_EV)) /
92,661✔
772
    (*E_electron + MASS_ELECTRON_EV);
92,661✔
773
  double rn = uniform_distribution(-1., 1., seed);
92,661✔
774
  *mu_electron = (rn + beta) / (rn * beta + 1.0);
92,661✔
775

776
  // Sample the scattering angle of the positron
777
  beta = std::sqrt(*E_positron * (*E_positron + 2.0 * MASS_ELECTRON_EV)) /
92,661✔
778
         (*E_positron + MASS_ELECTRON_EV);
92,661✔
779
  rn = uniform_distribution(-1., 1., seed);
92,661✔
780
  *mu_positron = (rn + beta) / (rn * beta + 1.0);
92,661✔
781
}
92,661✔
782

783
void PhotonInteraction::atomic_relaxation(int i_shell, Particle& p) const
25,432,633✔
784
{
785
  // Return if no atomic relaxation data is present or if the binding energy is
786
  // larger than the incident particle energy
787
  if (!has_atomic_relaxation_ || shells_[i_shell].binding_energy > p.E())
25,432,633!
UNCOV
788
    return;
×
789

790
  // Stack for unprocessed holes left by transitioning electrons
791
  int n_holes = 0;
25,432,633✔
792
  array<int, MAX_STACK_SIZE> holes;
25,432,633✔
793

794
  // Push the initial hole onto the stack
795
  holes[n_holes++] = i_shell;
25,432,633✔
796

797
  while (n_holes > 0) {
132,742,562✔
798
    // Pop the next hole off the stack
799
    int i_hole = holes[--n_holes];
107,309,929✔
800
    const auto& shell {shells_[i_hole]};
107,309,929✔
801

802
    // If no transitions, assume fluorescent photon from captured free electron
803
    if (shell.transitions.empty()) {
107,309,929✔
804
      Direction u = isotropic_direction(p.current_seed());
65,448,669✔
805
      double E = shell.binding_energy;
65,448,669✔
806
      p.create_secondary(p.wgt(), u, E, ParticleType::photon());
65,448,669✔
807
      continue;
65,448,669✔
808
    }
65,448,669✔
809

810
    // Sample transition
811
    double c = -prn(p.current_seed());
41,861,260✔
812
    int i_trans;
41,861,260✔
813
    for (i_trans = 0; i_trans < shell.transitions.size(); ++i_trans) {
1,005,900,933!
814
      c += shell.transitions[i_trans].probability;
1,005,900,933✔
815
      if (c > 0)
1,005,900,933✔
816
        break;
817
    }
818
    const auto& transition = shell.transitions[i_trans];
41,861,260✔
819

820
    // Sample angle isotropically
821
    Direction u = isotropic_direction(p.current_seed());
41,861,260✔
822

823
    // Push the hole created by the electron transitioning to the photoelectron
824
    // hole onto the stack
825
    holes[n_holes++] = transition.primary_subshell;
41,861,260✔
826

827
    if (transition.secondary_subshell != -1) {
41,861,260✔
828
      // Non-radiative transition -- Auger/Coster-Kronig effect
829

830
      // Push the hole left by emitted auger electron onto the stack
831
      holes[n_holes++] = transition.secondary_subshell;
40,016,036✔
832

833
      // Create auger electron
834
      p.create_secondary(
40,016,036✔
835
        p.wgt(), u, transition.energy, ParticleType::electron());
40,016,036✔
836
    } else {
837
      // Radiative transition -- get X-ray energy
838

839
      // Create fluorescent photon
840
      p.create_secondary(p.wgt(), u, transition.energy, ParticleType::photon());
1,845,224✔
841
    }
842
  }
843
}
844

845
//==============================================================================
846
// Non-member functions
847
//==============================================================================
848

849
std::pair<double, double> klein_nishina(double alpha, uint64_t* seed)
20,057,844✔
850
{
851
  double alpha_out, mu;
20,057,844✔
852
  double beta = 1.0 + 2.0 * alpha;
20,057,844✔
853
  if (alpha < 3.0) {
20,057,844✔
854
    // Kahn's rejection method
855
    double t = beta / (beta + 8.0);
19,380,253✔
856
    double x;
32,740,932✔
857
    while (true) {
32,740,932✔
858
      if (prn(seed) < t) {
32,740,932✔
859
        // Left branch of flow chart
860
        double r = uniform_distribution(0.0, 2.0, seed);
5,395,237✔
861
        x = 1.0 + alpha * r;
5,395,237✔
862
        if (prn(seed) < 4.0 / x * (1.0 - 1.0 / x)) {
5,395,237✔
863
          mu = 1 - r;
2,984,749✔
864
          break;
2,984,749✔
865
        }
866
      } else {
867
        // Right branch of flow chart
868
        x = beta / (1.0 + 2.0 * alpha * prn(seed));
27,345,695✔
869
        mu = 1.0 + (1.0 - x) / alpha;
27,345,695✔
870
        if (prn(seed) < 0.5 * (mu * mu + 1.0 / x))
27,345,695✔
871
          break;
872
      }
873
    }
874
    alpha_out = alpha / x;
19,380,253✔
875

876
  } else {
877
    // Koblinger's direct method
878
    double gamma = 1.0 - std::pow(beta, -2);
677,591✔
879
    double s =
677,591✔
880
      prn(seed) * (4.0 / alpha + 0.5 * gamma +
677,591✔
881
                    (1.0 - (1.0 + beta) / (alpha * alpha)) * std::log(beta));
677,591✔
882
    if (s <= 2.0 / alpha) {
677,591✔
883
      // For first term, x = 1 + 2ar
884
      // Therefore, a' = a/(1 + 2ar)
885
      alpha_out = alpha / (1.0 + 2.0 * alpha * prn(seed));
100,620✔
886
    } else if (s <= 4.0 / alpha) {
576,971✔
887
      // For third term, x = beta/(1 + 2ar)
888
      // Therefore, a' = a(1 + 2ar)/beta
889
      alpha_out = alpha * (1.0 + 2.0 * alpha * prn(seed)) / beta;
99,195✔
890
    } else if (s <= 4.0 / alpha + 0.5 * gamma) {
477,776✔
891
      // For fourth term, x = 1/sqrt(1 - gamma*r)
892
      // Therefore, a' = a*sqrt(1 - gamma*r)
893
      alpha_out = alpha * std::sqrt(1.0 - gamma * prn(seed));
123,786✔
894
    } else {
895
      // For third term, x = beta^r
896
      // Therefore, a' = a/beta^r
897
      alpha_out = alpha / std::pow(beta, prn(seed));
353,990✔
898
    }
899

900
    // Calculate cosine of scattering angle based on basic relation
901
    mu = 1.0 + 1.0 / alpha - 1.0 / alpha_out;
677,591✔
902
  }
903
  return {alpha_out, mu};
20,057,844✔
904
}
905

906
void free_memory_photon()
8,287✔
907
{
908
  data::elements.clear();
8,287✔
909
  data::compton_profile_pz.resize({0});
8,287✔
910
  data::ttb_e_grid.resize({0});
8,287✔
911
  data::ttb_k_grid.resize({0});
8,287✔
912
}
8,287✔
913

914
} // namespace openmc
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc