diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index dd2db069dde..c30d507398e 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -84,6 +84,9 @@ class Discrete : public Distribution { //! \param seed Pseudorandom number seed pointer //! \return Sampled value double sample(uint64_t* seed) const override; + //! Calculate the probability density function (PDF) at a given value + //! \param x The value at which to evaluate the PDF + //! \return The value of the PDF at the given point double get_pdf(double x) const; double integral() const override { return di_.integral(); }; diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 03a687b3875..74267064769 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -67,25 +67,98 @@ class FilterBinIter { //! \param p The particle being tracked void score_collision_tally(Particle& p); -//! Score tallies using the next event estimator. -// +//! Score tallies using the point estimator. //! This is triggered after every collision. -// //! \param p The particle being tracked void score_point_tally(Particle& p); +//! Retrieve the position and exclusion sphere (R0) of a detector. +//! +//! This function fetches the spatial position [x, y, z] and the exclusion sphere radius (R0) +//! of the point detector for a specified tally index. It ensures that the tally contains exactly +//! three spatial coordinates and an R0 value. +//! +//! \param det_pos Array to store the detector position [x, y, z] and exclusion sphere radius R0. +//! \param i_tally Index of the tally for which the detector position and R0 are required. void get_det_pos(double (&det_pos)[4], int i_tally); +//! Score to point tally using data from a source site. +//! +//! This function calculates and scores flux for active point tallies +//! based on the properties of a given source site. It accounts for the source's +//! spatial position, directional distribution, and type to compute the appropriate +//! contribution to the tally. The function handles various source angle distributions, +//! including isotropic and polar-azimuthal types. +//! +//! \param src Pointer to the source site containing particle and source information. +//! \throws RuntimeError If an unexpected angle distribution type is encountered. void score_point_tally_from_source(const SourceSite* src); +//! Calculate the total mean free path (MFP) for a ghost particle along a specified distance. +//! +//! This function computes the total mean free path (MFP) for a ghost particle traveling +//! from a collision point to a detector located at a given distance. The direction of travel +//! is towards the point detector. The cross-section used in the calculation is computed +//! based on the ghost particle's energy, which remains constant throughout the traversal. +//! The function accounts for shielding by iteratively advancing the particle to boundaries +//! and summing contributions to the MFP along the path. +//! +//! \param ghost_particle The particle whose MFP is being calculated, with its direction towards the point detector. +//! \param total_distance The distance [cm] from the collision point to the detector. +//! \return The total mean free path for the particle over the specified distance. +//! \throws RuntimeError If the particle encounters unexpected behavior during boundary crossing. double get_MFP(Particle& ghost_particle, double total_distance); +//! Calculate the probability density function (PDF) for elastic and inelastic scattering to a point detector. +//! +//! This function computes the PDF for a particle scattering (elastic or inelastic) to a detector +//! at a specified position. All computations are performed fully relativistically, including energy, +//! momentum, and transformations between the center-of-mass (CM) frame and the lab frame. The +//! Jacobian used in the calculations corresponds to a Lorentz transformation, ensuring accurate +//! mapping between frames of reference. +//! +//! Ghost particles are created to represent possible scattering trajectories. Each ghost particle +//! is initialized with the correct energy and direction as if the particle had scattered towards +//! the point detector. +//! +//! In the case of inelastic scattering, the post-collision energy of the outgoing neutron is set in the +//! CM frame, leading to an effective mass (\(m_4\)) for the residual nucleus. For inelastic collisions, +//! \(m_4 \neq m_2\), as \(m_4\) reflects the excitation state of the residual nucleus. + +//! \param det_pos The position of the point detector in [x, y, z, R0]. +//! \param p The incident particle being analyzed for scattering. +//! \param mu_cm A vector to store the cosine of the scattering angles in the CM frame. +//! \param Js A vector to store the Jacobians for transforming the distributions to the lab frame using Lorentz transformations. +//! \param ghost_particles A vector to store the ghost particles representing scattering trajectories. +//! Each ghost particle is initialized with the correct energy and direction +//! towards the point detector. +//! \param E3k_cm_given The center-of-mass kinetic energy of one outgoing particle, if specified. +//! If set to -1, it is calculated based on kinematics. void get_pdf_to_point_elastic(double det_pos[4], Particle& p, std::vector& mu_cm, std::vector& Js, std::vector& ghost_particles, double E3k_cm_given = -1); -void get_pdf_to_point_inelastic(double det_pos[4], Particle& p, - std::vector& mu_cm, std::vector& Js, - std::vector& ghost_particles, double E3_cm); -void get_mu_cm_inelastic(double det_pos[4], double E_out, Particle& p); +//! Score a ghost particle contribution to a point tally. +//! +//! This function calculates and scores the contribution of a ghost particle +//! to a specified point tally. The scoring process includes computing the flux +//! at the detector location. +//! It accounts for cases where the particle path is fully or partially within the exclusion +//! sphere (\(R_0\)) or outside it. The scoring respects particle type, ensuring only +//! neutrons are considered. +//! +//! \param ghost_p The ghost particle being scored. +//! \param pdf_lab The probability density function value in the lab frame for the particle direction. +//! \param i_tally The index of the tally to which the ghost particle contribution is scored. +//! +//! \note Only continuous energy mode is currently supported; multi-group mode is not implemented. void score_ghost_particle(Particle& ghost_p, double pdf_lab, int i_tally); -// +//! Perform a Lorentz boost of a four-vector. +//! +//! This function boosts a four-vector \( B \) (in the lab frame) into the rest frame +//! defined by another four-vector \( A \). The result of the boost is stored in \( X \), +//! representing \( B \) in the rest frame of \( A \). +//! +//! \param A The four-vector defining the rest frame [E, px, py, pz]. +//! \param B The four-vector to be boosted [E, px, py, pz]. +//! \param X The resulting four-vector after the Lorentz boost [E, px, py, pz] (rest frame). void boostf(double A[4], double B[4], double X[4]); + //! Analog tallies are triggered at every collision, not every event. // //! \param p The particle being tracked diff --git a/src/physics.cpp b/src/physics.cpp index 678ff8f67ac..165d4d7898c 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -1207,7 +1207,7 @@ void score_fission_neutron(int i_tally, int i_nuclide, const Reaction& rx, auto& ghost_p = ghost_particles[index]; double pdf_lab = pdfs_lab[index]; score_ghost_particle(ghost_p, pdf_lab, i_tally); - // calculate shielding + } // for loop on ghost particles } diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 0aff6bf0b13..4641f89f44f 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2728,7 +2728,7 @@ double get_MFP(Particle& ghost_particle, double total_distance) } remaining_distance -= advance_distance; ghost_particle.time() += - advance_distance / ghost_particle.speed(); // not reletevistic + advance_distance / ghost_particle.speed(); ghost_particle.event_cross_surface(); ghost_particle.event_calculate_xs(); ghost_particle.boundary() = distance_to_boundary(ghost_particle);