diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst
index 9f8eb84d80e..be9e10429da 100644
--- a/docs/source/methods/random_ray.rst
+++ b/docs/source/methods/random_ray.rst
@@ -994,6 +994,114 @@ random ray and Monte Carlo, however.
   develop the scattering source by way of inactive batches before beginning
   active batches.
 
+
+.. _usersguide_fixed_source_first_collision_source_methods:
+
+----------------------------------
+First Collision Source Method
+----------------------------------
+
+In cases where the fixed source is not a well-defined volumetric source (e.g. point source),
+there is the need to preprocess the source into volumetric sources to be computed by the
+random ray solver. One possible way is through the First Collision Source Method (FCS), which works as a 
+preconditioner of the source, distributing the source contribution throughout the domain. This is
+not a new method (`Alcouffe <Alcouffe-1989_>`_), and has been recently implemented in other 
+neutron transport codes (`Ragusa <Ragusa-2016_>`_, `Falabino <Falabino-2022_>`_).
+
+The FCS works by generating uncollided neutron angular fluxes :math:`\psi^{\text{un}}_{g}` at the source that travels through
+the domain interacting with the media, being naturally attenuated in the process. Neutrons that experience
+any collision are treated as first collided neutrons and will be used to estimate the volumetric neutron 
+fixed source at that cell. Once the FCS preprocess stage is complete, the fixed volumetric source will be
+added to each iteration of the random ray solver. The remaining uncollided flux that has not interacted 
+at any point of the preprocessing stage is added to the final solution at the end of the neutron transport code.
+
+The FCS has a very similar mathematical formulation to the regular Random Ray method. The formulation for
+the method with flat sources will be provided. The neutron transport equation for the uncollided rays can be
+described as in Equation :eq:`moc_final`, however without any pre-existing source terms:
+
+.. math::
+    :label: fcs_moc_final
+
+    \psi_g^{\text{un}}(s, \mathbf{\Omega}) = \psi_g^{\text{un}}(\mathbf{r_0}, \mathbf{\Omega}) e^{-\int_0^s ds^\prime \Sigma_{t_g}(s^\prime)}.
+
+The analytical solution for this ODE is a simple attenuation problem:
+
+.. math::
+    :label: fcs_moc_sol
+
+    \psi_g^{\text{un}}(s) = \psi_g(0) e^{-\Sigma_{t,i,g} s} .
+
+For convenience, we can also write this equation in terms of the incoming and
+outgoing angular flux (:math:`\psi_g^{in}` and :math:`\psi_g^{out}`) 
+
+.. math::
+    :label: fcs_fsr_attenuation_in_out
+
+    \psi_g^{out} = \psi_g^{in} e^{-\Sigma_{t,i,g} \ell_r}.
+
+
+Equation :eq:`fcs_fsr_attenuation_in_out` can be rearranged in terms of :math:`\Delta \psi_{r,g}`:
+
+.. math::
+    :label: fcs_difference
+    
+    \Delta\psi_{r,g}^{\text{un}} =\psi_{r,g}^{in} - \psi_{r,g}^{out} = \psi_{r,g}^{in}\big(1-e^{-\Sigma_{t,i,g}l_r}\big) .
+    
+The average angular flux can be computed substitute :eq:`average` into :eq:`fcs_moc_final`
+and obtain:
+
+.. math::
+    :label: fcs_average_solved
+
+    \overline{\psi}_{r,i,g}^{\text{un}} = - \frac{\psi_{r,g}^{out} - \psi_{r,g}^{in}}{\ell_r \Sigma_{t,i,g}}.
+
+Which can also be described in terms of :math:`\Delta \psi_{r,g}`:
+
+.. math::
+    :label: fcs_average_solved_difference
+
+    \overline{\psi}_{r,i,g}^{\text{un}} = \frac{\Delta \psi_{r,g}^{\text{un}}}{\ell_r \Sigma_{t,i,g}}.
+
+Similarly to Equation :eq:`discretized`, the scalar flux in cell :math:`i` can be defined as the summation
+of the contributions of the angular fluxes traveling through it. However, in this method, there is the need
+to differentiate how the volume is estimated. In the Random Ray method, the rays are generated uniformly
+across the domain and the volume estimation is unbiased. Nonetheless, angular fluxes have a specific and non-well-distributed
+birth location.
+
+.. math::
+    :label: fcs_discretized
+
+    \phi^{\text{un}}(i,g) = \frac{\int_{V_i} \int_{4\pi} \psi(r, \Omega) d\Omega d\mathbf{r}}{\int_{V_i} d\mathbf{r}} = \overline{\overline{\psi}}_{i,g} \approx \frac{\sum\limits_{r=1}^{N_i} \ell_r \overline{\psi}_{r,i,g}^{un}}{V_i} .
+
+To avoid bias in the volume estimation, it is added an initial volume estimation stage
+to compute the volume of each cell before running FCS, using the same method presented in 
+the Random Ray method. 
+
+In the case of flat sources, the volume estimation can be improved with the random ray solver, 
+allowing a fast and rough initial estimate. However, linear sources require a more careful procedure.
+Linear source requires calculating angular flux moments, that depend on the estimated
+cell's spatial moments and centroids. If a poor estimation is made, inaccurate gradients will be
+carried throughout the calculation and impact the accuracy of the method. 
+
+Substituting Equation :eq:`fcs_average_solved_difference` 
+into Equation :eq:`fcs_discretized`, we obtain the equation for the uncollided angular flux:
+
+.. math::
+    :label: fcs_uncollided_flux
+
+    \phi^{\text{un}}(i,g) =  \frac{\sum\limits_{r=1}^{N_i} \Delta \psi_{r,i,g}^{\text{un}}}{\Sigma_{t,i,g} V_i} .
+
+The fixed-source term can be calculated as the first collided flux that undergoes scattering events:
+
+.. math::
+    :label: fcs_first_collided_flux
+
+        Q_\text{fixed}(i,g) = \phi^{\text{FCS}}(i,g) = \sum_{g'}^{G}\Sigma_s(i,g,g') \phi^{\text{un}} (i,g') .
+
+The fixed source :math:`Q_\text{fixed}` will be treated as a fixed volumetric source for all remaining Random Ray
+iterations, as presented in Eq :eq:`fixed_source_update`.
+
+
 ---------------------------
 Fundamental Sources of Bias
 ---------------------------
@@ -1048,6 +1156,9 @@ in random ray particle transport are:
 .. _Cosgrove-2023: https://doi.org/10.1080/00295639.2023.2270618
 .. _Ferrer-2016: https://doi.org/10.13182/NSE15-6
 .. _Gunow-2018: https://dspace.mit.edu/handle/1721.1/119030
+.. _Alcouffe-1989: https://doi.org/10.13182/NSE90-A23749
+.. _Ragusa-2016: https://www.osti.gov/biblio/1364492
+.. _Falabino-2022: https://doi.org/10.1016/j.jcp.2022.111156
 
 .. only:: html
 
diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst
index 117d5e23fb5..50fec11f8f9 100644
--- a/docs/source/usersguide/random_ray.rst
+++ b/docs/source/usersguide/random_ray.rst
@@ -475,6 +475,41 @@ as::
 which will greatly improve the quality of the linear source term in 2D
 simulations.
 
+-----------------------------
+First Collision Source Method
+-----------------------------
+
+The First Collision Source method (FCS) is supported within the fixed source random ray
+solver. The preprocessing stage can be toggle by setting the ``first_collision_source``
+field in the :attr:`openmc.Settings.random_ray` dictionary to ``True`` as::
+
+    settings.random_ray['first_collision_source'] = True
+
+The method has two input variables: the number of rays for initial volume estimation
+and number of uncollided rays used in FCS.
+
+The user can define the amount of uncollided rays used in the FCS mode by setting the
+``first_collision_rays`` field in the :attr:`openmc.Settings.random_ray` dictionary 
+to ``value`` as::
+
+        settings.random_ray['first_collision_rays'] = 1000
+
+If not provided, the solver will run the default procedure to iteratively generate more
+rays until:
+
+* All cells are hit.
+* No extra cells were hit with additional rays.
+* Maximum value reached (default :math:`= 1e7` rays).
+
+The user can also define the amount of rays for initial volume estimation by setting the
+``first_collision_volume_rays`` field in the :attr:`openmc.Settings.random_ray` dictionary 
+to ``value`` as::
+    
+        settings.random_ray['first_collision_volume_rays'] = 2000
+
+If not provided, the solver will use the same amount of rays used for a regular batch in
+the random ray solver (:attr:`settings::n_particles`).
+
 ---------------------------------
 Fixed Source and Eigenvalue Modes
 ---------------------------------
diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
index e70c803de2f..797c3c18249 100644
--- a/include/openmc/distribution.h
+++ b/include/openmc/distribution.h
@@ -56,6 +56,7 @@ class DiscreteIndex {
   // Properties
   const vector<double>& prob() const { return prob_; }
   const vector<size_t>& alias() const { return alias_; }
+  const vector<double>& prob_actual() const { return prob_actual_; }
   double integral() const { return integral_; }
 
 private:
@@ -63,6 +64,7 @@ class DiscreteIndex {
                         //!< mapped to alias method table
   vector<size_t> alias_; //!< Alias table
   double integral_;      //!< Integral of distribution
+  vector<double> prob_actual_; //!< actual probability before the Vose algorithm
 
   //! Normalize distribution so that probabilities sum to unity
   void normalize();
@@ -91,6 +93,7 @@ class Discrete : public Distribution {
   const vector<double>& x() const { return x_; }
   const vector<double>& prob() const { return di_.prob(); }
   const vector<size_t>& alias() const { return di_.alias(); }
+  const vector<double>& prob_actual() const { return di_.prob_actual(); }
 
 private:
   vector<double> x_; //!< Possible outcomes
diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
index 164148cce10..78cb6047cb2 100644
--- a/include/openmc/particle_data.h
+++ b/include/openmc/particle_data.h
@@ -51,6 +51,7 @@ struct SourceSite {
   ParticleType particle;
   int64_t parent_id;
   int64_t progeny_id;
+  int source_id;
 };
 
 //! State of a particle used for particle track files
diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h
index 33c5661dcd4..3819810e7dd 100644
--- a/include/openmc/random_ray/flat_source_domain.h
+++ b/include/openmc/random_ray/flat_source_domain.h
@@ -113,6 +113,19 @@ class FlatSourceDomain {
   virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const;
   double compute_fixed_source_normalization_factor() const;
 
+  virtual void compute_first_collided_flux();
+  virtual void normalize_uncollided_scalar_flux(double number_of_particles);
+  virtual void update_volume_uncollided_flux();
+  virtual void compute_first_collided_external_source();
+  virtual void compute_uncollided_scalar_flux();
+  int64_t check_fsr_hits();
+  virtual void uncollided_moments();
+  virtual void batch_reset_fc();
+  virtual double evaluate_uncollided_flux_at_point(
+    Position r, int64_t sr, int g) const;
+  virtual double evaluate_external_source_at_point(
+    Position r, int64_t sr, int g) const;
+
   //----------------------------------------------------------------------------
   // Static Data members
   static bool volume_normalized_flux_tallies_;
@@ -144,6 +157,8 @@ class FlatSourceDomain {
   vector<float> scalar_flux_new_;
   vector<float> source_;
   vector<float> external_source_;
+  vector<float> scalar_uncollided_flux_;
+  vector<float> scalar_first_collided_flux_;
 
 protected:
   //----------------------------------------------------------------------------
diff --git a/include/openmc/random_ray/linear_source_domain.h b/include/openmc/random_ray/linear_source_domain.h
index 5010ffddd6f..eb18ca30a4e 100644
--- a/include/openmc/random_ray/linear_source_domain.h
+++ b/include/openmc/random_ray/linear_source_domain.h
@@ -41,6 +41,18 @@ class LinearSourceDomain : public FlatSourceDomain {
   void flux_swap() override;
   double evaluate_flux_at_point(Position r, int64_t sr, int g) const override;
 
+  void compute_first_collided_external_source() override;
+  void compute_uncollided_scalar_flux() override;
+  void compute_first_collided_flux() override;
+  void normalize_uncollided_scalar_flux(double number_of_particles) override;
+  void update_volume_uncollided_flux() override;
+  void uncollided_moments() override;
+  void batch_reset_fc() override;
+  double evaluate_uncollided_flux_at_point(
+    Position r, int64_t sr, int g) const override;
+  double evaluate_external_source_at_point(
+    Position r, int64_t sr, int g) const override;
+
   //----------------------------------------------------------------------------
   // Public Data members
 
@@ -48,6 +60,11 @@ class LinearSourceDomain : public FlatSourceDomain {
   vector<MomentArray> flux_moments_old_;
   vector<MomentArray> flux_moments_new_;
   vector<MomentArray> flux_moments_t_;
+  // First Collided Method
+  vector<MomentArray> flux_moments_uncollided_;
+  vector<MomentArray> flux_moments_first_collided_;
+  vector<MomentArray> external_source_gradients_;
+
   vector<Position> centroid_;
   vector<Position> centroid_iteration_;
   vector<Position> centroid_t_;
diff --git a/include/openmc/random_ray/random_ray.h b/include/openmc/random_ray/random_ray.h
index 913a9af4a75..c8a8f27b2fa 100644
--- a/include/openmc/random_ray/random_ray.h
+++ b/include/openmc/random_ray/random_ray.h
@@ -20,7 +20,7 @@ class RandomRay : public Particle {
   //----------------------------------------------------------------------------
   // Constructors
   RandomRay();
-  RandomRay(uint64_t ray_id, FlatSourceDomain* domain);
+  RandomRay(uint64_t ray_id, FlatSourceDomain* domain, bool uncollided_ray);
 
   //----------------------------------------------------------------------------
   // Methods
@@ -28,10 +28,10 @@ class RandomRay : public Particle {
   void attenuate_flux(double distance, bool is_active);
   void attenuate_flux_flat_source(double distance, bool is_active);
   void attenuate_flux_linear_source(double distance, bool is_active);
-
-  void initialize_ray(uint64_t ray_id, FlatSourceDomain* domain);
+  void event_advance_ray_first_collided();
+  void initialize_ray(uint64_t ray_id, FlatSourceDomain* domain,bool uncollided_ray);
   uint64_t transport_history_based_single_ray();
-
+  
   //----------------------------------------------------------------------------
   // Static data members
   static double distance_inactive_;          // Inactive (dead zone) ray length
@@ -39,9 +39,16 @@ class RandomRay : public Particle {
   static unique_ptr<Source> ray_source_;     // Starting source for ray sampling
   static RandomRaySourceShape source_shape_; // Flag for linear source
 
+  static bool first_collision_source_;
+  static int first_collision_rays_;
+  static int first_collision_volume_rays_;
+
+  static bool no_volume_calc; // Flag for FCS flux calculation 
+
   //----------------------------------------------------------------------------
   // Public data members
   vector<float> angular_flux_;
+  vector<float> angular_flux_initial_;
 
 private:
   //----------------------------------------------------------------------------
@@ -50,6 +57,7 @@ class RandomRay : public Particle {
   vector<MomentArray> delta_moments_;
 
   int negroups_;
+
   FlatSourceDomain* domain_ {nullptr}; // pointer to domain that has flat source
                                        // data needed for ray transport
   double distance_travelled_ {0};
diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h
index c1d47821d7a..52db67a3990 100644
--- a/include/openmc/random_ray/random_ray_simulation.h
+++ b/include/openmc/random_ray/random_ray_simulation.h
@@ -27,9 +27,13 @@ class RandomRaySimulation {
   void print_results_random_ray(uint64_t total_geometric_intersections,
     double avg_miss_rate, int negroups, int64_t n_source_regions,
     int64_t n_external_source_regions) const;
+  void first_collision_source_simulation();
 
   //----------------------------------------------------------------------------
   // Data members
+  // Initial volume estimation timer - First Collided Source Method  
+  double time_volume_fc;
+
 private:
   // Contains all flat source region data
   unique_ptr<FlatSourceDomain> domain_;
diff --git a/include/openmc/timer.h b/include/openmc/timer.h
index d928aad4560..5cffa82cc51 100644
--- a/include/openmc/timer.h
+++ b/include/openmc/timer.h
@@ -32,6 +32,7 @@ extern Timer time_event_surface_crossing;
 extern Timer time_event_collision;
 extern Timer time_event_death;
 extern Timer time_update_src;
+extern Timer time_first_collided;
 
 } // namespace simulation
 
diff --git a/openmc/__init__.py b/openmc/__init__.py
index 566d287068f..3c594fffa6a 100644
--- a/openmc/__init__.py
+++ b/openmc/__init__.py
@@ -40,5 +40,7 @@
 
 from . import examples
 
+#__version__ = '0.14.1-dev'
+
 
 __version__ = importlib.metadata.version("openmc")
diff --git a/openmc/lib/core.py b/openmc/lib/core.py
index a9a549fa05a..cfd3fdd8a58 100644
--- a/openmc/lib/core.py
+++ b/openmc/lib/core.py
@@ -26,7 +26,8 @@ class _SourceSite(Structure):
                 ('surf_id', c_int),
                 ('particle', c_int),
                 ('parent_id', c_int64),
-                ('progeny_id', c_int64)]
+                ('progeny_id', c_int64),
+                ('source_id', c_int)]
 
 
 # Define input type for numpy arrays that will be passed into C++ functions
diff --git a/openmc/settings.py b/openmc/settings.py
index ce97f138f24..1343ba4036c 100644
--- a/openmc/settings.py
+++ b/openmc/settings.py
@@ -164,6 +164,19 @@ class Settings:
             cm/cm^3. When disabled, flux tallies will be reported in units
             of cm (i.e., total distance traveled by neutrons in the spatial
             tally region).
+        :first_collision_source:
+            Boolean that indicates if first collided source method (FSC) is 
+            used to initialize the external source input. Default is 'False'.
+        :first_collision_rays:
+            Number of rays (int) used to generate the first collided source.
+            If not provided, the method will automatically run enough rays to
+            adequately converge the first collided source, via an iterative 
+            doubling process. If provided, the method will run the exact
+            prescribed amount.
+        :first_collision_volume_rays:
+            Number of rays (int) used to estimate the initial volume for the
+            first collied source calculation. If not provided, it will use the
+            same number of rays as the main random ray solver stage.
 
         .. versionadded:: 0.15.0
     resonance_scattering : dict
@@ -1096,6 +1109,14 @@ def random_ray(self, random_ray: dict):
                                ('flat', 'linear', 'linear_xy'))
             elif key == 'volume_normalized_flux_tallies':
                 cv.check_type('volume normalized flux tallies', random_ray[key], bool)
+            elif key == 'first_collision_source':
+                cv.check_type('first_collision_source', random_ray[key], bool)
+            elif key == 'first_collision_rays':
+                cv.check_type('first_collision_rays', random_ray[key], int)
+                cv.check_greater_than('first_collision_rays',random_ray[key], 0)
+            elif key == 'first_collision_volume_rays':
+                cv.check_type('first_collision_volume_rays', random_ray[key], int)
+                cv.check_greater_than('first_collision_volume_rays',random_ray[key], 0)                                
             else:
                 raise ValueError(f'Unable to set random ray to "{key}" which is '
                                  'unsupported by OpenMC')
@@ -1895,6 +1916,12 @@ def _random_ray_from_xml_element(self, root):
                     self.random_ray['volume_normalized_flux_tallies'] = (
                         child.text in ('true', '1')
                     )
+                elif child.tag == 'first_collision_source':
+                    self.random_ray['first_collision_source'] = (
+                        child.text in ('true', '1')
+                    )
+                elif child.tag in ('first_collision_rays', 'first_collision_volume_rays'):
+                    self.random_ray[child.tag] = int(child.text) 
 
     def to_xml_element(self, mesh_memo=None):
         """Create a 'settings' element to be written to an XML file.
diff --git a/src/distribution.cpp b/src/distribution.cpp
index 3026630b335..0d87fce801e 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -46,6 +46,9 @@ void DiscreteIndex::init_alias()
 {
   normalize();
 
+  // record user input normalized distribution prob_actual for Random Ray
+  prob_actual_ = prob_; 
+
   // The initialization and sampling method is based on Vose
   // (DOI: 10.1109/32.92917)
   // Vectors for large and small probabilities based on 1/n
diff --git a/src/initialize.cpp b/src/initialize.cpp
index cc1eac9cf35..e5de9b9ea39 100644
--- a/src/initialize.cpp
+++ b/src/initialize.cpp
@@ -155,7 +155,7 @@ void initialize_mpi(MPI_Comm intracomm)
 
   // Create bank datatype
   SourceSite b;
-  MPI_Aint disp[10];
+  MPI_Aint disp[11];
   MPI_Get_address(&b.r, &disp[0]);
   MPI_Get_address(&b.u, &disp[1]);
   MPI_Get_address(&b.E, &disp[2]);
@@ -166,14 +166,15 @@ void initialize_mpi(MPI_Comm intracomm)
   MPI_Get_address(&b.particle, &disp[7]);
   MPI_Get_address(&b.parent_id, &disp[8]);
   MPI_Get_address(&b.progeny_id, &disp[9]);
-  for (int i = 9; i >= 0; --i) {
+  MPI_Get_address(&b.source_id, &disp[10]);
+  for (int i = 10; i >= 0; --i) {
     disp[i] -= disp[0];
   }
 
-  int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1};
+  int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1};
   MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
-    MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG};
-  MPI_Type_create_struct(10, blocks, disp, types, &mpi::source_site);
+    MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG, MPI_INT};
+  MPI_Type_create_struct(11, blocks, disp, types, &mpi::source_site);
   MPI_Type_commit(&mpi::source_site);
 }
 #endif // OPENMC_MPI
diff --git a/src/output.cpp b/src/output.cpp
index 5fdbea1304e..2b1104aa929 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -32,6 +32,7 @@
 #include "openmc/nuclide.h"
 #include "openmc/plot.h"
 #include "openmc/random_ray/flat_source_domain.h"
+#include "openmc/random_ray/random_ray.h"
 #include "openmc/reaction.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
@@ -376,13 +377,16 @@ void print_build_info()
 
 void print_columns()
 {
-  if (settings::entropy_on) {
+  if (RandomRay::first_collision_source_){
+    fmt::print("  Batch       Rays       Total Source Regions Discovered\n"
+               "  ======      =======    =================================\n");
+  } else if (settings::entropy_on) {
     fmt::print("  Bat./Gen.      k       Entropy         Average k \n"
                "  =========   ========   ========   ====================\n");
   } else {
     fmt::print("  Bat./Gen.      k            Average k\n"
                "  =========   ========   ====================\n");
-  }
+  } 
 }
 
 //==============================================================================
diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp
index 8b6cda93f2a..57bc9062133 100644
--- a/src/random_ray/flat_source_domain.cpp
+++ b/src/random_ray/flat_source_domain.cpp
@@ -54,10 +54,14 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
   // Initialize element-wise arrays
   scalar_flux_new_.assign(n_source_elements_, 0.0);
   scalar_flux_final_.assign(n_source_elements_, 0.0);
-  source_.resize(n_source_elements_);
+  scalar_uncollided_flux_.assign(n_source_elements_, 0.0);
+  scalar_first_collided_flux_.assign(n_source_elements_, 0.0);
+  source_.assign(n_source_elements_, 0.0);
   external_source_.assign(n_source_elements_, 0.0);
   tally_task_.resize(n_source_elements_);
   volume_task_.resize(n_source_regions_);
+  //
+  // scalar_flux_new_per_source_.assign(n_source_elements_, 0.0);
 
   if (settings::run_mode == RunMode::EIGENVALUE) {
     // If in eigenvalue mode, set starting flux to guess of unity
@@ -110,8 +114,17 @@ void FlatSourceDomain::batch_reset()
   // Reset scalar fluxes, iteration volume tallies, and region hit flags to
   // zero
   parallel_fill<float>(scalar_flux_new_, 0.0f);
+  parallel_fill<int>(was_hit_, 0);
   parallel_fill<double>(volume_, 0.0);
+}
+
+void FlatSourceDomain::batch_reset_fc()
+{
+  parallel_fill<float>(scalar_flux_new_, 0.0f);
   parallel_fill<int>(was_hit_, 0);
+  if (RandomRay::no_volume_calc) {
+    parallel_fill<double>(volume_t_, 0.0);
+  }
 }
 
 void FlatSourceDomain::accumulate_iteration_flux()
@@ -122,6 +135,22 @@ void FlatSourceDomain::accumulate_iteration_flux()
   }
 }
 
+int64_t FlatSourceDomain::check_fsr_hits()
+{
+  int64_t n_hits = 0;
+#pragma omp parallel for reduction(+ : n_hits)
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+
+    // Check if this cell was hit this iteration
+    int was_cell_hit = was_hit_[sr];
+    if (was_cell_hit) {
+      n_hits++;
+    }
+  }
+
+  return n_hits;
+}
+
 // Compute new estimate of scattering + fission sources in each source region
 // based on the flux estimate from the previous iteration.
 void FlatSourceDomain::update_neutron_source(double k_eff)
@@ -200,7 +229,8 @@ void FlatSourceDomain::normalize_scalar_flux_and_volumes(
   double volume_normalization_factor =
     1.0 / (total_active_distance_per_iteration * simulation::current_batch);
 
-// Normalize scalar flux to total distance travelled by all rays this iteration
+  // Normalize scalar flux to total distance travelled by all rays this
+  // iteration
 #pragma omp parallel for
   for (int64_t e = 0; e < scalar_flux_new_.size(); e++) {
     scalar_flux_new_[e] *= normalization_factor;
@@ -215,6 +245,112 @@ void FlatSourceDomain::normalize_scalar_flux_and_volumes(
   }
 }
 
+void FlatSourceDomain::compute_uncollided_scalar_flux()
+{
+  // Temperature and angle indices, if using multiple temperature
+  const int t = 0;
+  const int a = 0;
+
+#pragma omp parallel for
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+    int was_cell_hit = was_hit_[sr];
+    int material = material_[sr];
+
+    for (int g = 0; g < negroups_; g++) {
+      int64_t idx = (sr * negroups_) + g;
+
+      if (was_cell_hit) {
+        float sigma_t = data::mg.macro_xs_[material].get_xs(
+          MgxsType::TOTAL, g, nullptr, nullptr, nullptr, t, a);
+        scalar_uncollided_flux_[idx] = scalar_flux_new_[idx] / (sigma_t);
+      }
+    }
+  }
+}
+
+// Normalize Uncollided Flux
+void FlatSourceDomain::normalize_uncollided_scalar_flux(
+  double number_of_particles)
+{
+  // multiply by simulation volume
+  float normalization_factor = (1.0) / number_of_particles;
+
+  // Determine Source_total Scailing factor if first collided
+  double user_external_source_strength = 0.0;
+  for (auto& ext_source : model::external_sources) {
+    user_external_source_strength += ext_source->strength();
+  }
+
+  float source_scailing_factor =
+    (user_external_source_strength * normalization_factor);
+
+#pragma omp parallel for
+  for (int64_t e = 0; e < scalar_uncollided_flux_.size(); e++) {
+    scalar_uncollided_flux_[e] *= source_scailing_factor;
+  }
+}
+
+// Compute First Collided flux
+void FlatSourceDomain::compute_first_collided_flux()
+{
+  const int t = 0;
+  const int a = 0;
+
+#pragma omp parallel for
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+    int material = material_[sr];
+
+    for (int e_out = 0; e_out < negroups_; e_out++) {
+      float sigma_t = data::mg.macro_xs_[material].get_xs(
+        MgxsType::TOTAL, e_out, nullptr, nullptr, nullptr, t, a);
+      float scatter_fixed_source = 0.0f;
+
+      for (int e_in = 0; e_in < negroups_; e_in++) {
+        float scalar_flux = scalar_uncollided_flux_[sr * negroups_ + e_in];
+
+        float sigma_s = data::mg.macro_xs_[material].get_xs(
+          MgxsType::NU_SCATTER, e_in, &e_out, nullptr, nullptr, t, a);
+        scatter_fixed_source += sigma_s * scalar_flux;
+      }
+
+      scalar_first_collided_flux_[sr * negroups_ + e_out] =
+        scatter_fixed_source / sigma_t;
+    }
+  }
+}
+
+void FlatSourceDomain::uncollided_moments()
+{ }
+
+// Normalize First Collided Flux by volume and assign as fixed source
+// compute_first_collided_external_source
+void FlatSourceDomain::compute_first_collided_external_source()
+{
+#pragma omp parallel for
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+    // check here
+    double volume = simulation_volume_ * volume_[sr];
+
+    if (volume > 0.0) {
+      for (int g = 0; g < negroups_; g++) {
+        external_source_[sr * negroups_ + g] =
+          (scalar_first_collided_flux_[sr * negroups_ + g] / (volume));
+      }
+    }
+  }
+}
+
+void FlatSourceDomain::update_volume_uncollided_flux()
+{
+#pragma omp parallel for
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+    double volume = volume_[sr] * simulation_volume_;
+    for (int g = 0; g < negroups_; g++) {
+      scalar_uncollided_flux_[sr * negroups_ + g] /= volume;
+    }
+  }
+}
+
 // Combine transport flux contributions and flat source contributions from the
 // previous iteration to generate this iteration's estimate of scalar flux.
 int64_t FlatSourceDomain::add_source_to_scalar_flux()
@@ -566,12 +702,17 @@ void FlatSourceDomain::random_ray_tally()
     // solves, but useful in fixed source solves for returning the flux shape
     // with a magnitude that makes sense relative to the fixed source strength.
     double volume = volume_[sr] * simulation_volume_;
-
     double material = material_[sr];
+    float flux = 0.0f;
+
     for (int g = 0; g < negroups_; g++) {
       int idx = sr * negroups_ + g;
-      double flux = scalar_flux_new_[idx] * source_normalization_factor;
-
+      if (RandomRay::first_collision_source_) {
+        flux =
+          (scalar_flux_new_[idx] + (scalar_uncollided_flux_[idx] / volume));
+      } else {
+        flux = (scalar_flux_new_[idx] * source_normalization_factor);
+      }
       // Determine numerical score value
       for (auto& task : tally_task_[idx]) {
         double score;
@@ -751,11 +892,39 @@ void FlatSourceDomain::all_reduce_replicated_source_regions()
 #endif
 }
 
+// avoid add scalar_uncollided_flux if not using fixed source - first collided
 double FlatSourceDomain::evaluate_flux_at_point(
   Position r, int64_t sr, int g) const
 {
-  return scalar_flux_final_[sr * negroups_ + g] /
-         (settings::n_batches - settings::n_inactive);
+  if (RandomRay::first_collision_source_) {
+    return (scalar_flux_final_[sr * negroups_ + g] /
+              (settings::n_batches - settings::n_inactive) +
+            scalar_uncollided_flux_[sr * negroups_ + g]);
+  } else {
+    // add uncollided flux if First Collided method is used
+    return (scalar_flux_final_[sr * negroups_ + g] /
+            (settings::n_batches - settings::n_inactive));
+  }
+}
+
+double FlatSourceDomain::evaluate_uncollided_flux_at_point(
+  Position r, int64_t sr, int g) const
+{
+  if (RandomRay::first_collision_source_) {
+    return scalar_uncollided_flux_[sr * negroups_ + g];
+  } else {
+    return 0;
+  }
+}
+
+double FlatSourceDomain::evaluate_external_source_at_point(
+  Position r, int64_t sr, int g) const
+{
+  if (RandomRay::first_collision_source_) {
+    return external_source_[sr * negroups_ + g];
+  } else {
+    return 0;
+  }
 }
 
 // Outputs all basic material, FSR ID, multigroup flux, and
@@ -906,6 +1075,37 @@ void FlatSourceDomain::output_to_vtk() const
       std::fwrite(&total_fission, sizeof(float), 1, plot);
     }
 
+    // Plot fixed source
+    for (int g = 0; g < negroups_; g++) {
+      std::fprintf(plot, "SCALARS fixed_source_group_%d float\n", g);
+      std::fprintf(plot, "LOOKUP_TABLE default\n");
+      for (int i = 0; i < Nx * Ny * Nz; i++) {
+        int64_t fsr = voxel_indices[i];
+        int mat = material_[fsr];
+        float sigma_t = data::mg.macro_xs_[mat].get_xs(
+          MgxsType::TOTAL, g, nullptr, nullptr, nullptr, 0, 0);
+        float f_source =
+          evaluate_external_source_at_point(voxel_positions[i], fsr, g);
+        f_source *= sigma_t;
+        f_source = convert_to_big_endian<float>(f_source);
+        std::fwrite(&f_source, sizeof(float), 1, plot);
+      }
+    }
+
+    // Plot multigroup uncollided flux data
+    for (int g = 0; g < negroups_; g++) {
+      std::fprintf(plot, "SCALARS Uncollided_flux_group_%d float\n", g);
+      std::fprintf(plot, "LOOKUP_TABLE default\n");
+      for (int i = 0; i < Nx * Ny * Nz; i++) {
+        int64_t fsr = voxel_indices[i];
+        int64_t source_element = fsr * negroups_ + g;
+        float uncollided_flux =
+          evaluate_uncollided_flux_at_point(voxel_positions[i], fsr, g);
+        uncollided_flux = convert_to_big_endian<float>(uncollided_flux);
+        std::fwrite(&uncollided_flux, sizeof(float), 1, plot);
+      }
+    }
+
     std::fclose(plot);
   }
 }
@@ -914,7 +1114,7 @@ void FlatSourceDomain::apply_external_source_to_source_region(
   Discrete* discrete, double strength_factor, int64_t source_region)
 {
   const auto& discrete_energies = discrete->x();
-  const auto& discrete_probs = discrete->prob();
+  const auto& discrete_probs = discrete->prob_actual();
 
   for (int e = 0; e < discrete_energies.size(); e++) {
     int g = data::mg.get_group_index(discrete_energies[e]);
diff --git a/src/random_ray/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp
index 1603ec24a4e..62539e2b908 100644
--- a/src/random_ray/linear_source_domain.cpp
+++ b/src/random_ray/linear_source_domain.cpp
@@ -28,6 +28,10 @@ LinearSourceDomain::LinearSourceDomain() : FlatSourceDomain()
   flux_moments_t_.assign(n_source_elements_, {0.0, 0.0, 0.0});
   // Source gradients given by M inverse multiplied by source moments
   source_gradients_.assign(n_source_elements_, {0.0, 0.0, 0.0});
+  external_source_gradients_.assign(n_source_elements_, {0.0, 0.0, 0.0});
+  // First Collided method
+  flux_moments_uncollided_.assign(n_source_elements_, {0.0, 0.0, 0.0});
+  flux_moments_first_collided_.assign(n_source_elements_, {0.0, 0.0, 0.0});
 
   centroid_.assign(n_source_regions_, {0.0, 0.0, 0.0});
   centroid_iteration_.assign(n_source_regions_, {0.0, 0.0, 0.0});
@@ -50,6 +54,43 @@ void LinearSourceDomain::batch_reset()
   }
 }
 
+void LinearSourceDomain::batch_reset_fc()
+{
+  FlatSourceDomain::batch_reset_fc();
+#pragma omp parallel for
+  for (int64_t se = 0; se < n_source_elements_; se++) {
+    flux_moments_new_[se] = {0.0, 0.0, 0.0};
+  }
+  if (RandomRay::no_volume_calc) {
+#pragma omp parallel for
+    for (int64_t sr = 0; sr < n_source_regions_; sr++) {
+      centroid_iteration_[sr] = {0.0, 0.0, 0.0};
+      centroid_[sr] = {0.0, 0.0, 0.0};
+      centroid_t_[sr] = {0.0, 0.0, 0.0};
+      mom_matrix_[sr] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+      mom_matrix_t_[sr] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+    }
+  }
+}
+
+void LinearSourceDomain::compute_first_collided_external_source()
+{
+  FlatSourceDomain::compute_first_collided_external_source();
+#pragma omp parallel for
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+
+    double volume = simulation_volume_ * volume_[sr];
+
+    // Normalize by volume and attribute it as external_source
+    if (volume > 0.0) {
+      for (int g = 0; g < negroups_; g++) {
+        external_source_gradients_[sr * negroups_ + g] =
+          (flux_moments_first_collided_[sr * negroups_ + g] * (1.0 / volume));
+      }
+    }
+  }
+}
+
 void LinearSourceDomain::update_neutron_source(double k_eff)
 {
   simulation::time_update_src.start();
@@ -103,7 +144,7 @@ void LinearSourceDomain::update_neutron_source(double k_eff)
         (scatter_flat + fission_flat * inverse_k_eff) / sigma_t;
 
       // Compute the linear source terms
-      if (simulation::current_batch > 2) {
+      if (simulation::current_batch > 10) {
         source_gradients_[sr * negroups_ + e_out] =
           invM * ((scatter_linear + fission_linear * inverse_k_eff) / sigma_t);
       }
@@ -115,6 +156,7 @@ void LinearSourceDomain::update_neutron_source(double k_eff)
 #pragma omp parallel for
     for (int se = 0; se < n_source_elements_; se++) {
       source_[se] += external_source_[se];
+      source_gradients_[se] += external_source_gradients_[se];
     }
   }
 
@@ -128,7 +170,7 @@ void LinearSourceDomain::normalize_scalar_flux_and_volumes(
   double volume_normalization_factor =
     1.0 / (total_active_distance_per_iteration * simulation::current_batch);
 
-// Normalize flux to total distance travelled by all rays this iteration
+  // Normalize flux to total distance travelled by all rays this iteration
 #pragma omp parallel for
   for (int64_t e = 0; e < scalar_flux_new_.size(); e++) {
     scalar_flux_new_[e] *= normalization_factor;
@@ -153,6 +195,15 @@ void LinearSourceDomain::normalize_scalar_flux_and_volumes(
   }
 }
 
+void LinearSourceDomain::compute_uncollided_scalar_flux()
+{
+#pragma omp parallel for
+  for (int64_t e = 0; e < scalar_flux_new_.size(); e++) {
+    scalar_uncollided_flux_[e] = scalar_flux_new_[e];
+    flux_moments_uncollided_[e] = flux_moments_new_[e];
+  }
+}
+
 int64_t LinearSourceDomain::add_source_to_scalar_flux()
 {
   int64_t n_hits = 0;
@@ -192,6 +243,7 @@ int64_t LinearSourceDomain::add_source_to_scalar_flux()
         // previous iteration, then we simply set the new scalar flux to be
         // equal to the contribution from the flat source alone.
         scalar_flux_new_[idx] = source_[idx];
+        // ??? flux_moments_new_[idx] = ????
       } else {
         // If the FSR was not hit this iteration, and it has never been hit in
         // any iteration (i.e., volume is zero), then we want to set this to 0
@@ -205,6 +257,88 @@ int64_t LinearSourceDomain::add_source_to_scalar_flux()
   return n_hits;
 }
 
+// Normalize Uncollided Flux
+void LinearSourceDomain::normalize_uncollided_scalar_flux(
+  double number_of_particles)
+{
+  // multiply by simulation volume
+  float normalization_factor = (1.0f) / number_of_particles;
+
+  // Determine Source_total Scailing factor if first collided
+  double user_external_source_strength = 0.0;
+  for (auto& ext_source : model::external_sources) {
+    user_external_source_strength += ext_source->strength();
+  }
+
+  double source_scailing_factor =
+    user_external_source_strength * normalization_factor;
+
+#pragma omp parallel for
+  for (int64_t e = 0; e < scalar_uncollided_flux_.size(); e++) {
+    scalar_uncollided_flux_[e] *= source_scailing_factor;
+    flux_moments_uncollided_[e] *= source_scailing_factor;
+  }
+}
+
+// Compute First Collided flux
+void LinearSourceDomain::compute_first_collided_flux()
+{
+  const int t = 0;
+  const int a = 0;
+
+#pragma omp parallel for
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+
+    int material = material_[sr];
+    MomentMatrix invM = mom_matrix_[sr].inverse();
+
+    for (int e_out = 0; e_out < negroups_; e_out++) {
+      float sigma_t = data::mg.macro_xs_[material].get_xs(
+        MgxsType::TOTAL, e_out, nullptr, nullptr, nullptr, t, a);
+
+      float scatter_flat = 0.0f;
+      MomentArray scatter_linear = {0.0, 0.0, 0.0};
+
+      for (int e_in = 0; e_in < negroups_; e_in++) {
+        // Handles for the flat and linear components of the flux
+        float flux_flat = scalar_uncollided_flux_[sr * negroups_ + e_in];
+        MomentArray flux_linear =
+          flux_moments_uncollided_[sr * negroups_ + e_in];
+
+        // Handles for cross sections
+        float sigma_s = data::mg.macro_xs_[material].get_xs(
+          MgxsType::NU_SCATTER, e_in, &e_out, nullptr, nullptr, t, a);
+
+        // Compute source terms for flat and linear components of the flux
+        scatter_flat += sigma_s * flux_flat;
+        scatter_linear += sigma_s * flux_linear;
+      }
+
+      // Compute the flat source term
+      scalar_first_collided_flux_[sr * negroups_ + e_out] =
+        (scatter_flat) / sigma_t;
+
+      // Compute the linear source terms
+      flux_moments_first_collided_[sr * negroups_ + e_out] =
+        invM * ((scatter_linear) / sigma_t);
+    }
+  }
+}
+
+void LinearSourceDomain::uncollided_moments()
+{
+
+#pragma omp parallel for
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+    MomentMatrix invM = mom_matrix_[sr].inverse();
+
+    for (int g = 0; g < negroups_; g++) {
+      int64_t idx = (sr * negroups_) + g;
+      flux_moments_uncollided_[idx] = invM * flux_moments_uncollided_[idx];
+    }
+  }
+}
+
 void LinearSourceDomain::flux_swap()
 {
   FlatSourceDomain::flux_swap();
@@ -232,9 +366,9 @@ void LinearSourceDomain::all_reduce_replicated_source_regions()
   // We are going to assume we can safely cast Position, MomentArray,
   // and MomentMatrix to contiguous arrays of doubles for the MPI
   // allreduce operation. This is a safe assumption as typically
-  // compilers will at most pad to 8 byte boundaries. If a new FP32 MomentArray
-  // type is introduced, then there will likely be padding, in which case this
-  // function will need to become more complex.
+  // compilers will at most pad to 8 byte boundaries. If a new FP32
+  // MomentArray type is introduced, then there will likely be padding, in
+  // which case this function will need to become more complex.
   if (sizeof(MomentArray) != 3 * sizeof(double) ||
       sizeof(MomentMatrix) != 6 * sizeof(double)) {
     fatal_error("Unexpected buffer padding in linear source domain reduction.");
@@ -254,6 +388,7 @@ void LinearSourceDomain::all_reduce_replicated_source_regions()
 double LinearSourceDomain::evaluate_flux_at_point(
   Position r, int64_t sr, int g) const
 {
+  // RR neutron flux
   float phi_flat = FlatSourceDomain::evaluate_flux_at_point(r, sr, g);
 
   Position local_r = r - centroid_[sr];
@@ -263,7 +398,49 @@ double LinearSourceDomain::evaluate_flux_at_point(
   MomentMatrix invM = mom_matrix_[sr].inverse();
   MomentArray phi_solved = invM * phi_linear;
 
+  if (RandomRay::first_collision_source_) {
+    phi_solved += flux_moments_uncollided_[sr * negroups_ + g];
+  }
+
   return phi_flat + phi_solved.dot(local_r);
 }
 
+double LinearSourceDomain::evaluate_uncollided_flux_at_point(
+  Position r, int64_t sr, int g) const
+{
+  // Uncollided flux 
+  float phi_flat = FlatSourceDomain::evaluate_uncollided_flux_at_point(r, sr, g);
+
+  Position local_r = r - centroid_[sr];
+  MomentArray phi_linear = flux_moments_uncollided_[sr * negroups_ + g];
+
+  return phi_flat + phi_linear.dot(local_r);
+}
+
+double LinearSourceDomain::evaluate_external_source_at_point(
+  Position r, int64_t sr, int g) const
+{
+  // External source 
+  float phi_flat = FlatSourceDomain::evaluate_external_source_at_point(r, sr, g);
+
+  Position local_r = r - centroid_[sr];
+  MomentArray phi_linear = external_source_gradients_[sr * negroups_ + g];
+
+  return phi_flat + phi_linear.dot(local_r);
+}
+
+void LinearSourceDomain::update_volume_uncollided_flux()
+{
+#pragma omp parallel for
+  for (int sr = 0; sr < n_source_regions_; sr++) {
+    double volume = volume_[sr] * simulation_volume_;
+    if (volume > 0.0) {
+      for (int g = 0; g < negroups_; g++) {
+        scalar_uncollided_flux_[sr * negroups_ + g] *= (1.0 / volume);
+        flux_moments_uncollided_[sr * negroups_ + g] *= (1.0 / volume);
+      }
+    }
+  }
+}
+
 } // namespace openmc
diff --git a/src/random_ray/moment_matrix.cpp b/src/random_ray/moment_matrix.cpp
index 0324a14943b..547d24b27c3 100644
--- a/src/random_ray/moment_matrix.cpp
+++ b/src/random_ray/moment_matrix.cpp
@@ -25,7 +25,7 @@ MomentMatrix MomentMatrix::inverse() const
 
   // Check if the determinant is zero
   double det = determinant();
-  if (det < std::abs(1.0e-10)) {
+  if (det < std::abs(1.0e-6)) {
     // Set the inverse to zero. In effect, this will
     // result in all the linear terms of the source becoming
     // zero, leaving just the flat source.
diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp
index a5bf6ec1060..7d4b050695e 100644
--- a/src/random_ray/random_ray.cpp
+++ b/src/random_ray/random_ray.cpp
@@ -183,11 +183,16 @@ double RandomRay::distance_inactive_;
 double RandomRay::distance_active_;
 unique_ptr<Source> RandomRay::ray_source_;
 RandomRaySourceShape RandomRay::source_shape_ {RandomRaySourceShape::FLAT};
+bool RandomRay::no_volume_calc = {false};
+bool RandomRay::first_collision_source_ {false};
+int RandomRay::first_collision_rays_ {-1};
+int RandomRay::first_collision_volume_rays_ {-1};
 
 RandomRay::RandomRay()
   : angular_flux_(data::mg.num_energy_groups_),
     delta_psi_(data::mg.num_energy_groups_),
-    negroups_(data::mg.num_energy_groups_)
+    negroups_(data::mg.num_energy_groups_),
+    angular_flux_initial_(data::mg.num_energy_groups_)
 {
   if (source_shape_ == RandomRaySourceShape::LINEAR ||
       source_shape_ == RandomRaySourceShape::LINEAR_XY) {
@@ -195,9 +200,11 @@ RandomRay::RandomRay()
   }
 }
 
-RandomRay::RandomRay(uint64_t ray_id, FlatSourceDomain* domain) : RandomRay()
+RandomRay::RandomRay(
+  uint64_t ray_id, FlatSourceDomain* domain, bool uncollided_ray)
+  : RandomRay()
 {
-  initialize_ray(ray_id, domain);
+  initialize_ray(ray_id, domain, uncollided_ray);
 }
 
 // Transports ray until termination criteria are met
@@ -205,7 +212,11 @@ uint64_t RandomRay::transport_history_based_single_ray()
 {
   using namespace openmc;
   while (alive()) {
-    event_advance_ray();
+    if (no_volume_calc) {
+      event_advance_ray_first_collided();
+    } else {
+      event_advance_ray();
+    }
     if (!alive())
       break;
     event_cross_surface();
@@ -272,6 +283,51 @@ void RandomRay::event_advance_ray()
   }
 }
 
+// Transports uncollided ray across a single region.
+void RandomRay::event_advance_ray_first_collided()
+{
+  // Find the distance to the nearest boundary
+  boundary() = distance_to_boundary(*this);
+  double distance = boundary().distance;
+
+  if (distance <= 0.0) {
+    mark_as_lost("Negative transport distance detected for particle " +
+                 std::to_string(id()));
+    return;
+  }
+  // For Uncollided/First Collided Flux, it is calculated the attenuation
+  // as the ray advance through the region and a check if the outcoming
+  // flux reaches the defined threshold.
+  distance_travelled_ += distance;
+  attenuate_flux(distance, true);
+
+  // Advance particle
+  for (int j = 0; j < n_coord(); ++j) {
+    coord(j).r += distance * coord(j).u;
+  }
+
+  // Terminate ray threshold
+  float ray_threshold = 1e-20f;
+
+  // Check if angular flux is not zero or above threshold
+  // Do not terminate the ray if any of these conditions are met
+  // terminate otherwise
+  bool angular_flux_below_threshold = true;
+  for (int g = 0; g < negroups_; g++) {
+    if (angular_flux_initial_[g] > 0) {
+      // calculate the attenuation of ray (kills if ratio below threshold)
+      float ratio = angular_flux_[g] / angular_flux_initial_[g];
+      if (ratio >= ray_threshold) {
+        angular_flux_below_threshold = false;
+        break;
+      }
+    }
+  }
+  if (angular_flux_below_threshold) {
+    wgt() = 0.0;
+  }
+}
+
 void RandomRay::attenuate_flux(double distance, bool is_active)
 {
   switch (source_shape_) {
@@ -347,6 +403,11 @@ void RandomRay::attenuate_flux_flat_source(double distance, bool is_active)
     for (int g = 0; g < negroups_; g++) {
       domain_->scalar_flux_new_[source_element + g] += delta_psi_[g];
     }
+    // Accomulate volume (ray distance) into this iteration's estimate
+    // of the source region's volume
+    if (!no_volume_calc) {
+      domain_->volume_[source_region] += distance;
+    }
 
     // If the source region hasn't been hit yet this iteration,
     // indicate that it now has
@@ -354,10 +415,6 @@ void RandomRay::attenuate_flux_flat_source(double distance, bool is_active)
       domain_->was_hit_[source_region] = 1;
     }
 
-    // Accomulate volume (ray distance) into this iteration's estimate
-    // of the source region's volume
-    domain_->volume_[source_region] += distance;
-
     // Tally valid position inside the source region (e.g., midpoint of
     // the ray) if not done already
     if (!domain_->position_recorded_[source_region]) {
@@ -365,7 +422,6 @@ void RandomRay::attenuate_flux_flat_source(double distance, bool is_active)
       domain_->position_[source_region] = midpoint;
       domain_->position_recorded_[source_region] = 1;
     }
-
     // Release lock
     domain_->lock_[source_region].unlock();
   }
@@ -441,6 +497,7 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
     // Compute linear source terms, spatial and directional (dir),
     // calculated from the source gradients dot product with local centroid
     // and direction, respectively.
+
     float spatial_source =
       domain_->source_[source_element + g] +
       rm_local.dot(domain->source_gradients_[source_element + g]);
@@ -491,6 +548,7 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
 
     // Accumulate deltas into the new estimate of source region flux for this
     // iteration
+
     for (int g = 0; g < negroups_; g++) {
       domain_->scalar_flux_new_[source_element + g] += delta_psi_[g];
       domain->flux_moments_new_[source_element + g] += delta_moments_[g];
@@ -500,10 +558,12 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
     // momement estimates into the running totals for the iteration for this
     // source region. The centroid and spatial momements estimates are scaled by
     // the ray segment length as part of length averaging of the estimates.
-    domain_->volume_[source_region] += distance;
-    domain->centroid_iteration_[source_region] += midpoint * distance;
-    moment_matrix_estimate *= distance;
-    domain->mom_matrix_[source_region] += moment_matrix_estimate;
+    if (!no_volume_calc) {
+      domain_->volume_[source_region] += distance;
+      domain->centroid_iteration_[source_region] += midpoint * distance;
+      moment_matrix_estimate *= distance;
+      domain->mom_matrix_[source_region] += moment_matrix_estimate;
+    }
 
     // If the source region hasn't been hit yet this iteration,
     // indicate that it now has
@@ -523,7 +583,9 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
   }
 }
 
-void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
+// updated
+void RandomRay::initialize_ray(
+  uint64_t ray_id, FlatSourceDomain* domain, bool uncollided_ray)
 {
   domain_ = domain;
 
@@ -543,12 +605,17 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
   init_particle_seeds(particle_seed, seeds());
   stream() = STREAM_TRACKING;
 
-  // Sample from ray source distribution
-  SourceSite site {ray_source_->sample(current_seed())};
+  // Sample ray from input source distribution
+  SourceSite site;
+  if (uncollided_ray) {
+    site = sample_external_source(current_seed());
+  } else {
+    site = ray_source_->sample(current_seed());
+  }
   site.E = lower_bound_index(
     data::mg.rev_energy_bins_.begin(), data::mg.rev_energy_bins_.end(), site.E);
   site.E = negroups_ - site.E - 1.;
-  this->from_source(&site);
+  from_source(&site);
 
   // Locate ray
   if (lowest_coord().cell == C_NONE) {
@@ -562,14 +629,34 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
       cell_born() = lowest_coord().cell;
   }
 
-  // Initialize ray's starting angular flux to starting location's isotropic
-  // source
-  int i_cell = lowest_coord().cell;
-  int64_t source_region_idx =
-    domain_->source_region_offsets_[i_cell] + cell_instance();
+  // Initialize angular flux spectrum from source input
+  if (uncollided_ray) {
+    Source* source_handle = model::external_sources[site.source_id].get();
+    IndependentSource* is = dynamic_cast<IndependentSource*>(source_handle);
+    Discrete* energy_external_source = dynamic_cast<Discrete*>(is->energy());
 
-  for (int g = 0; g < negroups_; g++) {
-    angular_flux_[g] = domain_->source_[source_region_idx * negroups_ + g];
+    const auto& discrete_energies = energy_external_source->x();
+    const auto& discrete_probs = energy_external_source->prob_actual();
+
+    std::fill(angular_flux_.begin(), angular_flux_.end(), 0.0f);
+    std::fill(angular_flux_initial_.begin(), angular_flux_initial_.end(), 0.0f);
+
+    for (int e = 0; e < discrete_energies.size(); e++) {
+      int g = data::mg.get_group_index(discrete_energies[e]);
+      angular_flux_[g] = discrete_probs[e];
+      angular_flux_initial_[g] = angular_flux_[g];
+    }
+  } else {
+    // Initialize ray's starting angular flux to starting location's isotropic
+    // source
+    int i_cell = lowest_coord().cell;
+    int64_t source_region_idx =
+      domain_->source_region_offsets_[i_cell] + cell_instance();
+
+    for (int g = 0; g < negroups_; g++) {
+      angular_flux_[g] = domain_->source_[source_region_idx * negroups_ + g];
+      // check for zero angular_flux_[g] causing ray to die?
+    }
   }
 }
 
diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp
index 4bc77645bcd..009809b434a 100644
--- a/src/random_ray/random_ray_simulation.cpp
+++ b/src/random_ray/random_ray_simulation.cpp
@@ -166,21 +166,22 @@ void validate_random_ray_inputs()
           "random ray mode");
       }
 
-      // Check for isotropic source
-      UnitSphereDistribution* angle_dist = is->angle();
-      Isotropic* id = dynamic_cast<Isotropic*>(angle_dist);
-      if (!id) {
-        fatal_error(
-          "Invalid source definition -- only isotropic external sources are "
-          "allowed in random ray mode.");
-      }
-
-      // Validate that a domain ID was specified
-      if (is->domain_ids().size() == 0) {
-        fatal_error("Fixed sources must be specified by domain "
-                    "id (cell, material, or universe) in random ray mode.");
+      if (!RandomRay::first_collision_source_) {
+        // Check for isotropic source
+        UnitSphereDistribution* angle_dist = is->angle();
+        Isotropic* id = dynamic_cast<Isotropic*>(angle_dist);
+        if (!id) {
+          fatal_error(
+            "Invalid source definition -- only isotropic external sources are "
+            "allowed in random ray mode.");
+        }
+
+        // Validate that a domain ID was specified
+        if (is->domain_ids().size() == 0) {
+          fatal_error("Fixed sources must be specified by domain "
+                      "id (cell, material, or universe) in random ray mode.");
+        }
       }
-
       // Check that a discrete energy distribution was used
       Distribution* d = is->energy();
       Discrete* dd = dynamic_cast<Discrete*>(d);
@@ -248,6 +249,8 @@ RandomRaySimulation::RandomRaySimulation()
     domain_ = make_unique<FlatSourceDomain>();
     break;
   case RandomRaySourceShape::LINEAR:
+    domain_ = make_unique<LinearSourceDomain>();
+    break;
   case RandomRaySourceShape::LINEAR_XY:
     domain_ = make_unique<LinearSourceDomain>();
     break;
@@ -259,11 +262,14 @@ RandomRaySimulation::RandomRaySimulation()
 void RandomRaySimulation::simulate()
 {
   if (settings::run_mode == RunMode::FIXED_SOURCE) {
-    // Transfer external source user inputs onto random ray source regions
-    domain_->convert_external_sources();
-    domain_->count_external_source_regions();
+    if (!RandomRay::first_collision_source_) {
+      // Transfer external source user inputs onto random ray source regions
+      domain_->convert_external_sources();
+      domain_->count_external_source_regions();
+    } else {
+      first_collision_source_simulation();
+    }
   }
-
   // Random ray power iteration loop
   while (simulation::current_batch < settings::n_batches) {
 
@@ -274,6 +280,11 @@ void RandomRaySimulation::simulate()
     // Reset total starting particle weight used for normalizing tallies
     simulation::total_weight = 1.0;
 
+    // Update external source if FIRST_COLLIDED_METHOD is used
+    if (RandomRay::first_collision_source_) {
+      domain_->compute_first_collided_external_source();
+    }
+
     // Update source term (scattering + fission)
     domain_->update_neutron_source(k_eff_);
 
@@ -288,7 +299,7 @@ void RandomRaySimulation::simulate()
 #pragma omp parallel for schedule(dynamic)                                     \
   reduction(+ : total_geometric_intersections_)
     for (int i = 0; i < simulation::work_per_rank; i++) {
-      RandomRay ray(i, domain_.get());
+      RandomRay ray(i, domain_.get(), false);
       total_geometric_intersections_ +=
         ray.transport_history_based_single_ray();
     }
@@ -338,6 +349,10 @@ void RandomRaySimulation::simulate()
     finalize_generation();
     finalize_batch();
   } // End random ray power iteration loop
+
+  if (RandomRay::first_collision_source_) {
+    domain_->update_volume_uncollided_flux();
+  }
 }
 
 void RandomRaySimulation::reduce_simulation_statistics()
@@ -416,8 +431,10 @@ void RandomRaySimulation::print_results_random_ray(
     fmt::print(
       " Total Iterations                  = {}\n", settings::n_batches);
     fmt::print(" Flat Source Regions (FSRs)        = {}\n", n_source_regions);
-    fmt::print(
-      " FSRs Containing External Sources  = {}\n", n_external_source_regions);
+    if (!RandomRay::first_collision_source_) {
+      fmt::print(
+        " FSRs Containing External Sources  = {}\n", n_external_source_regions);
+    }
     fmt::print(" Total Geometric Intersections     = {:.4e}\n",
       static_cast<double>(total_geometric_intersections));
     fmt::print("   Avg per Iteration               = {:.4e}\n",
@@ -435,6 +452,10 @@ void RandomRaySimulation::print_results_random_ray(
     header("Timing Statistics", 4);
     show_time("Total time for initialization", time_initialize.elapsed());
     show_time("Reading cross sections", time_read_xs.elapsed(), 1);
+    if (RandomRay::first_collision_source_) {
+      show_time("Volume estimation time", RandomRaySimulation::time_volume_fc);
+      show_time("First Collided Source time", time_first_collided.elapsed());
+    }
     show_time("Total simulation time", time_total.elapsed());
     show_time("Transport sweep only", time_transport.elapsed(), 1);
     show_time("Source update only", time_update_src.elapsed(), 1);
@@ -457,4 +478,121 @@ void RandomRaySimulation::print_results_random_ray(
   }
 }
 
+void RandomRaySimulation::first_collision_source_simulation()
+{
+  if (RandomRay::source_shape_ == RandomRaySourceShape::FLAT) {
+    header("FIRST COLLISION SOURCE METHOD - Flat source", 3);
+  } else {
+    header("FIRST COLLISION SOURCE METHOD - Linear source", 3);
+  }
+
+  simulation::time_first_collided.start();
+  simulation::current_batch = 1;
+  fmt::print("Initial volume estimation...");
+  // Volume estimation is necessary to scale the first collided source
+  // accordingly. Estimation of linear moments in linear source has direct
+  // impact into final solution accuracy.
+
+  // Check if number of volume rays were provided or set default
+  if (RandomRay::first_collision_volume_rays_ == -1) {
+    RandomRay::first_collision_volume_rays_ = settings::n_particles;
+  }
+
+  // Ray tracing - volume calculation
+#pragma omp parallel for
+  for (int i = 0; i < RandomRay::first_collision_volume_rays_; i++) {
+    RandomRay ray(i, domain_.get(), false);
+    ray.transport_history_based_single_ray();
+  }
+
+  // Normalizing volumes
+  domain_->normalize_scalar_flux_and_volumes(
+    RandomRay::first_collision_volume_rays_ * RandomRay::distance_active_);
+
+  // Reset parameters for First Collided Source Method
+  domain_->batch_reset_fc();
+  RandomRay::no_volume_calc = true;
+
+  // print volume estimation time
+  time_volume_fc = simulation::time_first_collided.elapsed();
+  fmt::print("completed.\n\n");
+
+  //==============================================================================
+  // First Collided Source - Transport Loop
+  // It will operate until meet any criteria (1-3):
+  // (1) There isn't new FSR hits from batch_first_collided (n-1) to the next
+  // (n) (2) Reached pre-set maximum n_uncollided_rays (3) Hit 100% of the FSRs
+
+  // Prepare column header
+  print_columns();
+
+  // Default values
+  int new_n_rays = 1000;
+  int n_rays_max = 10000000;
+
+  // initialize loop variables
+  int64_t n_hits_old = 0;
+  int64_t n_hits_new = 0;
+  int batch_first_collided = 1;
+  int old_n_rays = 0;
+
+  // Check if number of rays were provided and update limit for loop
+  if (RandomRay::first_collision_rays_ != -1) {
+    new_n_rays = RandomRay::first_collision_rays_;
+    n_rays_max = RandomRay::first_collision_rays_;
+  }
+
+  while (old_n_rays < n_rays_max) { // Condition (2)
+
+    // Ray tracing and attenuation
+#pragma omp parallel for
+    for (int i = old_n_rays; i < new_n_rays; i++) {
+      RandomRay ray(i, domain_.get(), true);
+      ray.transport_history_based_single_ray();
+    }
+
+    int64_t n_hits_new = domain_->check_fsr_hits();
+
+    // Print results
+    fmt::print(
+      "  {:6}   {:10}    {:}\n", batch_first_collided, new_n_rays, n_hits_new);
+
+    // Update values for next batch
+    batch_first_collided++;
+    old_n_rays = new_n_rays;
+    new_n_rays *= 2;
+
+    // BREAK STATEMENT
+    if (n_hits_new == n_hits_old) {
+      break; // Condition (1)
+    } else if (static_cast<double>(n_hits_new) >= domain_->n_source_regions_) {
+      break; // Condition (3)
+    }
+
+    n_hits_old = n_hits_new;
+  }
+
+  // Compute scalar_uncollided_flux
+  domain_->compute_uncollided_scalar_flux();
+
+  // Normalize scalar_uncollided_flux
+  domain_->normalize_uncollided_scalar_flux(old_n_rays);
+
+  // compute first_collided_fixed_source
+  domain_->compute_first_collided_flux();
+
+  // apply invM*uncollided flux
+  domain_->uncollided_moments();
+
+  // reset values for RandomRay iteration
+  domain_
+    ->batch_reset_fc(); // clean-up of key variables (preserves just volume_)
+  RandomRay::no_volume_calc = {false};
+  simulation::current_batch = 0; // garantee the first batch will be 1 in RR
+
+  // compute First Collided Method simulation time
+  simulation::time_first_collided.stop();
+  fmt::print("\n");
+}
+
 } // namespace openmc
diff --git a/src/settings.cpp b/src/settings.cpp
index db956abf5ca..06a31cb6b6b 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -25,6 +25,7 @@
 #include "openmc/plot.h"
 #include "openmc/random_lcg.h"
 #include "openmc/random_ray/random_ray.h"
+#include "openmc/random_ray/random_ray_simulation.h"
 #include "openmc/simulation.h"
 #include "openmc/source.h"
 #include "openmc/string_utils.h"
@@ -286,6 +287,26 @@ void get_run_parameters(pugi::xml_node node_base)
       FlatSourceDomain::volume_normalized_flux_tallies_ =
         get_node_value_bool(random_ray_node, "volume_normalized_flux_tallies");
     }
+    if (check_for_node(random_ray_node, "first_collision_source")) {
+      RandomRay::first_collision_source_ =
+        get_node_value_bool(random_ray_node, "first_collision_source");
+
+      if (check_for_node(random_ray_node, "first_collision_rays")) {
+        RandomRay::first_collision_rays_ =
+          std::stoi(get_node_value(random_ray_node, "first_collision_rays"));
+        if (RandomRay::first_collision_rays_ <= 0) {
+          fatal_error("Number of first collided rays must be greater than 0");
+        }
+      }
+      if (check_for_node(random_ray_node, "first_collision_volume_rays")) {
+        RandomRay::first_collision_volume_rays_ = std::stoi(
+          get_node_value(random_ray_node, "first_collision_volume_rays"));
+        if (RandomRay::first_collision_volume_rays_ <= 0) {
+          fatal_error("Number of rays for first collided initial volume "
+                      "estimation must be greater than 0");
+        }
+      }
+    }
   }
 }
 
diff --git a/src/source.cpp b/src/source.cpp
index 15fe8433ba5..1157a1cfc44 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -627,7 +627,7 @@ SourceSite sample_external_source(uint64_t* seed)
 
   // Sample source site from i-th source distribution
   SourceSite site {model::external_sources[i]->sample_with_constraints(seed)};
-
+  site.source_id = i;
   // If running in MG, convert site.E to group
   if (!settings::run_CE) {
     site.E = lower_bound_index(data::mg.rev_energy_bins_.begin(),
diff --git a/src/timer.cpp b/src/timer.cpp
index 6d692d4fbf6..7694f660d2e 100644
--- a/src/timer.cpp
+++ b/src/timer.cpp
@@ -27,6 +27,7 @@ Timer time_event_surface_crossing;
 Timer time_event_collision;
 Timer time_event_death;
 Timer time_update_src;
+Timer time_first_collided;
 
 } // namespace simulation
 
@@ -87,6 +88,7 @@ void reset_timers()
   simulation::time_event_collision.reset();
   simulation::time_event_death.reset();
   simulation::time_update_src.reset();
+  simulation::time_first_collided.reset();
 }
 
 } // namespace openmc
diff --git a/tests/regression_tests/random_ray_first_collision/__init__.py b/tests/regression_tests/random_ray_first_collision/__init__.py
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/tests/regression_tests/random_ray_first_collision/flat/inputs_true.dat b/tests/regression_tests/random_ray_first_collision/flat/inputs_true.dat
new file mode 100644
index 00000000000..0ab42bac23c
--- /dev/null
+++ b/tests/regression_tests/random_ray_first_collision/flat/inputs_true.dat
@@ -0,0 +1,245 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="source">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="source"/>
+    </material>
+    <material id="2" name="void">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="void"/>
+    </material>
+    <material id="3" name="absorber">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="absorber"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="infinite source region" universe="1"/>
+    <cell id="2" material="2" name="infinite void region" universe="2"/>
+    <cell id="3" material="3" name="infinite absorber region" universe="3"/>
+    <cell fill="4" id="4" universe="5"/>
+    <cell fill="5" id="5" name="full domain" region="1 -2 3 -4 5 -6" universe="6"/>
+    <lattice id="4">
+      <pitch>2.5 2.5 2.5</pitch>
+      <dimension>12 12 12</dimension>
+      <lower_left>0.0 0.0 0.0</lower_left>
+      <universes>
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 </universes>
+    </lattice>
+    <surface boundary="reflective" coeffs="0.0" id="1" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="2" type="x-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="3" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="4" type="y-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="5" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="6" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>90</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source particle="neutron" strength="3.14" type="independent">
+      <space type="box">
+        <parameters>0.0 0.0 0.0 5.0 5.0 5.0</parameters>
+      </space>
+      <energy type="discrete">
+        <parameters>100.0 1.0</parameters>
+      </energy>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <distance_active>500.0</distance_active>
+      <distance_inactive>100.0</distance_inactive>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
+        </space>
+      </source>
+      <volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
+      <source_shape>flat</source_shape>
+      <first_collision_source>True</first_collision_source>
+    </random_ray>
+  </settings>
+  <tallies>
+    <filter id="3" type="material">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="material">
+      <bins>2</bins>
+    </filter>
+    <filter id="1" type="material">
+      <bins>3</bins>
+    </filter>
+    <tally id="3" name="Source Tally">
+      <filters>3</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+    <tally id="2" name="Void Tally">
+      <filters>2</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+    <tally id="1" name="Absorber Tally">
+      <filters>1</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+  </tallies>
+</model>
diff --git a/tests/regression_tests/random_ray_first_collision/flat/results_true.dat b/tests/regression_tests/random_ray_first_collision/flat/results_true.dat
new file mode 100644
index 00000000000..4a51b64a825
--- /dev/null
+++ b/tests/regression_tests/random_ray_first_collision/flat/results_true.dat
@@ -0,0 +1,9 @@
+tally 1:
+5.237693E-01
+5.603795E-02
+tally 2:
+2.626249E-02
+1.379769E-04
+tally 3:
+1.832107E-03
+6.713236E-07
diff --git a/tests/regression_tests/random_ray_first_collision/linear/inputs_true.dat b/tests/regression_tests/random_ray_first_collision/linear/inputs_true.dat
new file mode 100644
index 00000000000..b3dff653b3e
--- /dev/null
+++ b/tests/regression_tests/random_ray_first_collision/linear/inputs_true.dat
@@ -0,0 +1,245 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="source">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="source"/>
+    </material>
+    <material id="2" name="void">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="void"/>
+    </material>
+    <material id="3" name="absorber">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="absorber"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="infinite source region" universe="1"/>
+    <cell id="2" material="2" name="infinite void region" universe="2"/>
+    <cell id="3" material="3" name="infinite absorber region" universe="3"/>
+    <cell fill="4" id="4" universe="5"/>
+    <cell fill="5" id="5" name="full domain" region="1 -2 3 -4 5 -6" universe="6"/>
+    <lattice id="4">
+      <pitch>2.5 2.5 2.5</pitch>
+      <dimension>12 12 12</dimension>
+      <lower_left>0.0 0.0 0.0</lower_left>
+      <universes>
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 </universes>
+    </lattice>
+    <surface boundary="reflective" coeffs="0.0" id="1" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="2" type="x-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="3" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="4" type="y-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="5" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="6" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>90</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source particle="neutron" strength="3.14" type="independent">
+      <space type="box">
+        <parameters>0.0 0.0 0.0 5.0 5.0 5.0</parameters>
+      </space>
+      <energy type="discrete">
+        <parameters>100.0 1.0</parameters>
+      </energy>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <distance_active>500.0</distance_active>
+      <distance_inactive>100.0</distance_inactive>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
+        </space>
+      </source>
+      <volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
+      <source_shape>linear</source_shape>
+      <first_collision_source>True</first_collision_source>
+    </random_ray>
+  </settings>
+  <tallies>
+    <filter id="3" type="material">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="material">
+      <bins>2</bins>
+    </filter>
+    <filter id="1" type="material">
+      <bins>3</bins>
+    </filter>
+    <tally id="3" name="Source Tally">
+      <filters>3</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+    <tally id="2" name="Void Tally">
+      <filters>2</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+    <tally id="1" name="Absorber Tally">
+      <filters>1</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+  </tallies>
+</model>
diff --git a/tests/regression_tests/random_ray_first_collision/linear/results_true.dat b/tests/regression_tests/random_ray_first_collision/linear/results_true.dat
new file mode 100644
index 00000000000..58587cf74fe
--- /dev/null
+++ b/tests/regression_tests/random_ray_first_collision/linear/results_true.dat
@@ -0,0 +1,9 @@
+tally 1:
+1.684539E-01
+9.991864E-01
+tally 2:
+-3.623011E-02
+3.710201E-04
+tally 3:
+2.184070E-03
+9.589758E-07
diff --git a/tests/regression_tests/random_ray_first_collision/test.py b/tests/regression_tests/random_ray_first_collision/test.py
new file mode 100644
index 00000000000..9481e83048c
--- /dev/null
+++ b/tests/regression_tests/random_ray_first_collision/test.py
@@ -0,0 +1,36 @@
+import os
+
+import numpy as np
+import openmc
+from openmc.utility_funcs import change_directory
+from openmc.examples import random_ray_three_region_cube
+import pytest
+
+from tests.testing_harness import TolerantPyAPITestHarness
+
+
+class MGXSTestHarness(TolerantPyAPITestHarness):
+    def _cleanup(self):
+        super()._cleanup()
+        f = 'mgxs.h5'
+        if os.path.exists(f):
+            os.remove(f)
+
+
+@pytest.mark.parametrize("shape", ["flat", "linear"])
+def test_random_ray_first_collision(shape):
+    with change_directory(shape):
+        openmc.reset_auto_ids()
+        model = random_ray_three_region_cube()
+        model.settings.random_ray['source_shape'] = shape
+        model.settings.random_ray['first_collision_source'] = True
+        strengths = [1.0] # Good - fast group appears largest (besides most thermal)
+        midpoints = [100.0]
+        energy_dist = openmc.stats.Discrete(x=midpoints,p=strengths)
+        lower_left_src = [0.0, 0.0, 0.0]
+        upper_right_src = [5.0, 5.0, 5.0]
+        spatial_distribution = openmc.stats.Box(lower_left_src, upper_right_src, only_fissionable=False)
+        source = openmc.IndependentSource(space=spatial_distribution, energy=energy_dist, strength = 3.14)
+        model.settings.source = [source]
+        harness = MGXSTestHarness('statepoint.10.h5', model)
+        harness.main()
\ No newline at end of file
diff --git a/tests/regression_tests/random_ray_first_collision_rays/__init__.py b/tests/regression_tests/random_ray_first_collision_rays/__init__.py
new file mode 100644
index 00000000000..e69de29bb2d
diff --git a/tests/regression_tests/random_ray_first_collision_rays/inputs_true.dat b/tests/regression_tests/random_ray_first_collision_rays/inputs_true.dat
new file mode 100644
index 00000000000..53d98f04f22
--- /dev/null
+++ b/tests/regression_tests/random_ray_first_collision_rays/inputs_true.dat
@@ -0,0 +1,246 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material id="1" name="source">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="source"/>
+    </material>
+    <material id="2" name="void">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="void"/>
+    </material>
+    <material id="3" name="absorber">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="absorber"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="infinite source region" universe="1"/>
+    <cell id="2" material="2" name="infinite void region" universe="2"/>
+    <cell id="3" material="3" name="infinite absorber region" universe="3"/>
+    <cell fill="4" id="4" universe="5"/>
+    <cell fill="5" id="5" name="full domain" region="1 -2 3 -4 5 -6" universe="6"/>
+    <lattice id="4">
+      <pitch>2.5 2.5 2.5</pitch>
+      <dimension>12 12 12</dimension>
+      <lower_left>0.0 0.0 0.0</lower_left>
+      <universes>
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+1 1 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+2 2 2 2 2 2 2 2 2 2 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 
+3 3 3 3 3 3 3 3 3 3 3 3 </universes>
+    </lattice>
+    <surface boundary="reflective" coeffs="0.0" id="1" type="x-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="2" type="x-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="3" type="y-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="4" type="y-plane"/>
+    <surface boundary="reflective" coeffs="0.0" id="5" type="z-plane"/>
+    <surface boundary="vacuum" coeffs="30.0" id="6" type="z-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>90</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source particle="neutron" strength="3.14" type="independent">
+      <space type="box">
+        <parameters>0.0 0.0 0.0 5.0 5.0 5.0</parameters>
+      </space>
+      <energy type="discrete">
+        <parameters>100.0 1.0</parameters>
+      </energy>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <distance_active>500.0</distance_active>
+      <distance_inactive>100.0</distance_inactive>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
+        </space>
+      </source>
+      <volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
+      <first_collision_source>True</first_collision_source>
+      <first_collision_rays>1246</first_collision_rays>
+      <first_collision_volume_rays>2000</first_collision_volume_rays>
+    </random_ray>
+  </settings>
+  <tallies>
+    <filter id="3" type="material">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="material">
+      <bins>2</bins>
+    </filter>
+    <filter id="1" type="material">
+      <bins>3</bins>
+    </filter>
+    <tally id="3" name="Source Tally">
+      <filters>3</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+    <tally id="2" name="Void Tally">
+      <filters>2</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+    <tally id="1" name="Absorber Tally">
+      <filters>1</filters>
+      <scores>flux</scores>
+      <estimator>tracklength</estimator>
+    </tally>
+  </tallies>
+</model>
diff --git a/tests/regression_tests/random_ray_first_collision_rays/results_true.dat b/tests/regression_tests/random_ray_first_collision_rays/results_true.dat
new file mode 100644
index 00000000000..481d3fa224d
--- /dev/null
+++ b/tests/regression_tests/random_ray_first_collision_rays/results_true.dat
@@ -0,0 +1,9 @@
+tally 1:
+5.179652E-01
+5.485091E-02
+tally 2:
+2.634285E-02
+1.388239E-04
+tally 3:
+1.832152E-03
+6.713564E-07
diff --git a/tests/regression_tests/random_ray_first_collision_rays/test.py b/tests/regression_tests/random_ray_first_collision_rays/test.py
new file mode 100644
index 00000000000..6bdfdaf5e61
--- /dev/null
+++ b/tests/regression_tests/random_ray_first_collision_rays/test.py
@@ -0,0 +1,35 @@
+import os
+
+import numpy as np
+import openmc
+from openmc.utility_funcs import change_directory
+from openmc.examples import random_ray_three_region_cube
+import pytest
+
+from tests.testing_harness import TolerantPyAPITestHarness
+
+
+class MGXSTestHarness(TolerantPyAPITestHarness):
+    def _cleanup(self):
+        super()._cleanup()
+        f = 'mgxs.h5'
+        if os.path.exists(f):
+            os.remove(f)
+
+
+def test_random_ray_first_collision_rays():
+    openmc.reset_auto_ids()
+    model = random_ray_three_region_cube()
+    model.settings.random_ray['first_collision_source'] = True
+    model.settings.random_ray['first_collision_rays'] = 1246
+    model.settings.random_ray['first_collision_volume_rays'] = 2000
+    strengths = [1.0] 
+    midpoints = [100.0]
+    energy_dist = openmc.stats.Discrete(x=midpoints,p=strengths)
+    lower_left_src = [0.0, 0.0, 0.0]
+    upper_right_src = [5.0, 5.0, 5.0]
+    spatial_distribution = openmc.stats.Box(lower_left_src, upper_right_src, only_fissionable=False)
+    source = openmc.IndependentSource(space=spatial_distribution, energy=energy_dist, strength = 3.14)
+    model.settings.source = [source]
+    harness = MGXSTestHarness('statepoint.10.h5', model)
+    harness.main()
\ No newline at end of file