diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9fe133a22e3..9e7aecefa49 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -355,6 +355,7 @@ endif()
 #===============================================================================
 
 list(APPEND libopenmc_SOURCES
+  src/angle_energy.cpp
   src/atomic_mass.cpp
   src/bank.cpp
   src/boundary_condition.cpp
@@ -451,6 +452,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_musurface.cpp
   src/tallies/filter_parent_nuclide.cpp
   src/tallies/filter_particle.cpp
+  src/tallies/filter_point.cpp  
   src/tallies/filter_particle_production.cpp
   src/tallies/filter_polar.cpp
   src/tallies/filter_reaction.cpp
@@ -461,6 +463,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_universe.cpp
   src/tallies/filter_weight.cpp
   src/tallies/filter_zernike.cpp
+  src/tallies/next_event_scoring.cpp
   src/tallies/tally.cpp
   src/tallies/tally_scoring.cpp
   src/tallies/trigger.cpp
diff --git a/include/openmc/angle_energy.h b/include/openmc/angle_energy.h
index 55deb5d415e..1ac8f6feb8e 100644
--- a/include/openmc/angle_energy.h
+++ b/include/openmc/angle_energy.h
@@ -27,12 +27,20 @@ class AngleEnergy {
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  virtual double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const = 0;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  virtual double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const = 0;
   virtual ~AngleEnergy() = default;
 };
 
+double get_jac_and_transform(
+  double E_in, double& mu, double& E_out, uint64_t* seed, double awr);
+
+double get_jac_and_transform_impl(
+  double E_com, double& mu, double& E_out, uint64_t* seed, double awr);
+
 } // namespace openmc
 
 #endif // OPENMC_ANGLE_ENERGY_H
diff --git a/include/openmc/chain.h b/include/openmc/chain.h
index e01f2a01c18..99a8749d54a 100644
--- a/include/openmc/chain.h
+++ b/include/openmc/chain.h
@@ -76,9 +76,11 @@ class DecayPhotonAngleEnergy : public AngleEnergy {
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   const Distribution* photon_energy_;
diff --git a/include/openmc/constants.h b/include/openmc/constants.h
index 0b425a673dc..85d6e1ebe93 100644
--- a/include/openmc/constants.h
+++ b/include/openmc/constants.h
@@ -296,9 +296,9 @@ enum class MgxsType {
 
 enum class TallyResult { VALUE, SUM, SUM_SQ, SUM_THIRD, SUM_FOURTH };
 
-enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT };
+enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT, POINT };
 
-enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION };
+enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION, NEXT_EVENT };
 
 enum class TallyEvent { SURFACE, LATTICE, KILL, SCATTER, ABSORB };
 
diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
index 5b632bbce53..52aba848a10 100644
--- a/include/openmc/particle_data.h
+++ b/include/openmc/particle_data.h
@@ -484,7 +484,7 @@ class GeometryState {
  *   Algorithms.” Annals of Nuclear Energy 113 (March 2018): 506–18.
  *   https://doi.org/10.1016/j.anucene.2017.11.032.
  */
-class ParticleData : public GeometryState {
+class ParticleData : virtual public GeometryState {
 private:
   //==========================================================================
   // Data members -- see public: below for descriptions
diff --git a/include/openmc/physics.h b/include/openmc/physics.h
index 2472d979939..9cc6d95dc02 100644
--- a/include/openmc/physics.h
+++ b/include/openmc/physics.h
@@ -85,7 +85,7 @@ void sample_fission_neutron(
 
 //! handles all reactions with a single secondary neutron (other than fission),
 //! i.e. level scattering, (n,np), (n,na), etc.
-void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p);
+void inelastic_scatter(int i_nuclide, const Reaction& rx, Particle& p);
 
 void sample_secondary_photons(Particle& p, int i_nuclide);
 
diff --git a/include/openmc/position.h b/include/openmc/position.h
index 5d291d26b95..92346e27ad4 100644
--- a/include/openmc/position.h
+++ b/include/openmc/position.h
@@ -204,6 +204,17 @@ inline Position operator/(double a, Position b)
   return b /= a;
 }
 
+inline bool operator<(Position a, Position b)
+{
+  if (a[0] != b[0])
+    return (a[0] < b[0]);
+  if (a[1] != b[1])
+    return (a[1] < b[1]);
+  if (a[2] != b[2])
+    return (a[2] < b[2]);
+  return false;
+}
+
 inline Position Position::reflect(Position n) const
 {
   const double projection = n.dot(*this);
diff --git a/include/openmc/ray.h b/include/openmc/ray.h
index 62e86b0d90f..b357ae8bbc4 100644
--- a/include/openmc/ray.h
+++ b/include/openmc/ray.h
@@ -1,14 +1,14 @@
 #ifndef OPENMC_RAY_H
 #define OPENMC_RAY_H
 
-#include "openmc/particle_data.h"
+#include "openmc/particle.h"
 #include "openmc/position.h"
 
 namespace openmc {
 
 // Base class that implements ray tracing logic, not necessarily through
 // defined regions of the geometry but also outside of it.
-class Ray : public GeometryState {
+class Ray : virtual public GeometryState {
 
 public:
   // Initialize from location and direction
@@ -24,7 +24,7 @@ class Ray : public GeometryState {
    * Traces the ray through the geometry, calling on_intersection
    * at every surface boundary.
    */
-  void trace();
+  void trace(double max_distance = INFTY);
 
   // Stops the ray and exits tracing when called from on_intersection
   void stop() { stop_ = true; }
@@ -32,6 +32,8 @@ class Ray : public GeometryState {
   // Sets the dist_ variable
   void compute_distance();
 
+  virtual void update_distance();
+
 protected:
   // Records how far the ray has traveled
   double traversal_distance_ {0.0};
@@ -46,5 +48,32 @@ class Ray : public GeometryState {
   unsigned event_counter_ {0};
 };
 
+class ParticleRay : public Ray, public Particle {
+
+public:
+  ParticleRay(
+    Position r, Direction u, ParticleType type_, double time_, double E_)
+    : Ray(r, u)
+  {
+    type() = type_;
+    time() = time_;
+    E() = E_;
+    E_last() = E_;
+    r_last() = r;
+  }
+
+  void on_intersection() override;
+
+  // Sets the dist_ variable
+  void update_distance() override;
+
+  const double& traversal_distance() const { return traversal_distance_; }
+  const double& traversal_mfp() const { return traversal_mfp_; }
+
+protected:
+  // Records how much mean free paths the ray traveled
+  double traversal_mfp_ {0.0};
+};
+
 } // namespace openmc
 #endif // OPENMC_RAY_H
diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h
index 79a6160e25b..ba464633650 100644
--- a/include/openmc/reaction_product.h
+++ b/include/openmc/reaction_product.h
@@ -60,9 +60,11 @@ class ReactionProduct {
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const;
 
   ParticleType particle_;      //!< Particle type
   EmissionMode emission_mode_; //!< Emission mode
diff --git a/include/openmc/secondary_correlated.h b/include/openmc/secondary_correlated.h
index 69b22981aff..df5bd9e9c35 100644
--- a/include/openmc/secondary_correlated.h
+++ b/include/openmc/secondary_correlated.h
@@ -53,9 +53,11 @@ class CorrelatedAngleEnergy : public AngleEnergy {
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
   // energy property
   vector<double>& energy() { return energy_; }
diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h
index b25352be9f0..22e55285edb 100644
--- a/include/openmc/secondary_kalbach.h
+++ b/include/openmc/secondary_kalbach.h
@@ -46,9 +46,11 @@ class KalbachMann : public AngleEnergy {
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   //! Outgoing energy/angle at a single incoming energy
diff --git a/include/openmc/secondary_nbody.h b/include/openmc/secondary_nbody.h
index 9d033a6b869..eb1b53b9dba 100644
--- a/include/openmc/secondary_nbody.h
+++ b/include/openmc/secondary_nbody.h
@@ -39,9 +39,11 @@ class NBodyPhaseSpace : public AngleEnergy {
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   int n_bodies_;      //!< Number of particles distributed
diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h
index 45d2c5260f1..54d791d8498 100644
--- a/include/openmc/secondary_thermal.h
+++ b/include/openmc/secondary_thermal.h
@@ -39,9 +39,11 @@ class CoherentElasticAE : public AngleEnergy {
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
@@ -74,9 +76,11 @@ class IncoherentElasticAE : public AngleEnergy {
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   double debye_waller_;
@@ -108,9 +112,11 @@ class IncoherentElasticAEDiscrete : public AngleEnergy {
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   const vector<double>& energy_;  //!< Energies at which cosines are tabulated
@@ -149,9 +155,11 @@ class IncoherentInelasticAEDiscrete : public AngleEnergy {
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   const vector<double>& energy_; //!< Incident energies
@@ -196,9 +204,11 @@ class IncoherentInelasticAE : public AngleEnergy {
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   //! Secondary energy/angle distribution
@@ -246,9 +256,11 @@ class MixedElasticAE : public AngleEnergy {
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
 private:
   CoherentElasticAE coherent_dist_;         //!< Coherent distribution
diff --git a/include/openmc/secondary_uncorrelated.h b/include/openmc/secondary_uncorrelated.h
index f895ae77f63..94301e7e16d 100644
--- a/include/openmc/secondary_uncorrelated.h
+++ b/include/openmc/secondary_uncorrelated.h
@@ -37,9 +37,11 @@ class UncorrelatedAngleEnergy : public AngleEnergy {
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
-  double sample_energy_and_pdf(
-    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
+  double sample_energy_and_pdf(double E_in, double mu, double& E_out,
+    uint64_t* seed, bool is_com, double awr) const override;
 
   // Accessors
   AngleDistribution& angle() { return angle_; }
diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h
index 9a6cf1b2131..3484ab707d4 100644
--- a/include/openmc/simulation.h
+++ b/include/openmc/simulation.h
@@ -36,12 +36,14 @@ extern "C" double
 extern double log_spacing; //!< lethargy spacing for energy grid searches
 extern "C" int n_lost_particles;   //!< cumulative number of lost particles
 extern "C" bool need_depletion_rx; //!< need to calculate depletion rx?
-extern "C" int restart_batch;      //!< batch at which a restart job resumed
-extern "C" bool satisfy_triggers;  //!< have tally triggers been satisfied?
-extern int ssw_current_file;       //!< current surface source file
-extern "C" int total_gen;          //!< total number of generations simulated
-extern double total_weight;        //!< Total source weight in a batch
-extern int64_t work_per_rank;      //!< number of particles per MPI rank
+extern bool
+  nonvacuum_boundary_present;     //!< Does the geometry contain non-vacuum b.c.
+extern "C" int restart_batch;     //!< batch at which a restart job resumed
+extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied?
+extern int ssw_current_file;      //!< current surface source file
+extern "C" int total_gen;         //!< total number of generations simulated
+extern double total_weight;       //!< Total source weight in a batch
+extern int64_t work_per_rank;     //!< number of particles per MPI rank
 
 extern const RegularMesh* entropy_mesh;
 extern const RegularMesh* ufs_mesh;
diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h
index 77b0d9f420d..34553f5dc1d 100644
--- a/include/openmc/tallies/filter.h
+++ b/include/openmc/tallies/filter.h
@@ -40,6 +40,7 @@ enum class FilterType {
   PARENT_NUCLIDE,
   PARTICLE,
   PARTICLE_PRODUCTION,
+  POINT,
   POLAR,
   REACTION,
   SPHERICAL_HARMONICS,
diff --git a/include/openmc/tallies/filter_point.h b/include/openmc/tallies/filter_point.h
new file mode 100644
index 00000000000..909f0baa27d
--- /dev/null
+++ b/include/openmc/tallies/filter_point.h
@@ -0,0 +1,55 @@
+#ifndef OPENMC_TALLIES_FILTER_POINT_H
+#define OPENMC_TALLIES_FILTER_POINT_H
+
+#include "openmc/position.h"
+#include "openmc/span.h"
+#include "openmc/tallies/filter.h"
+#include "openmc/vector.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins tally by point detectors
+//==============================================================================
+
+class PointFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~PointFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "point"; }
+  FilterType type() const override { return FilterType::POINT; }
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator,
+    FilterMatch& match) const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  const vector<std::pair<Position, double>>& detectors() const
+  {
+    return detectors_;
+  }
+
+  void set_detectors(span<std::pair<Position, double>> detectors);
+
+private:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  vector<std::pair<Position, double>> detectors_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_PARTICLE_H
diff --git a/include/openmc/tallies/next_event_scoring.h b/include/openmc/tallies/next_event_scoring.h
new file mode 100644
index 00000000000..59247d9b302
--- /dev/null
+++ b/include/openmc/tallies/next_event_scoring.h
@@ -0,0 +1,60 @@
+#ifndef OPENMC_TALLIES_NEXT_EVENT_SCORING_H
+#define OPENMC_TALLIES_NEXT_EVENT_SCORING_H
+
+#include "openmc/nuclide.h"
+#include "openmc/particle.h"
+#include "openmc/random_lcg.h"
+#include "openmc/ray.h"
+#include "openmc/tallies/filter.h"
+#include "openmc/tallies/tally.h"
+#include "openmc/tallies/tally_scoring.h"
+#include "openmc/thermal.h"
+
+namespace openmc {
+
+void score_point_tally_elastic(
+  Particle& p, int i_nuclide, const Reaction& rx, int i_product, Direction v_t);
+
+void score_point_tally_inelastic(
+  Particle& p, int i_nuclide, const Reaction& rx, int i_product, double yield);
+
+void score_point_tally_fission(
+  Particle& p, int i_nuclide, const Reaction& rx, int i_product);
+
+void score_point_tally_sab(Particle& p, int i_nuclide, const ThermalData& sab,
+  const NuclideMicroXS& micro);
+
+void score_point_tally_source(SourceSite& site, int source_index);
+
+template<typename PDF>
+void score_point_tally_impl(
+  const Position r, const ParticleType type, const double time, PDF pdffunc)
+{
+  for (auto& det : model::active_point_detectors) {
+    auto u = (det - r);
+    double total_distance = u.norm();
+    u /= total_distance;
+    double E;
+    double pdf = pdffunc(u, E);
+    if (pdf == 0.0)
+      continue;
+    auto p = ParticleRay(r, u, type, time, E);
+    p.Ray::trace(total_distance);
+    double distance = p.traversal_distance();
+    if (distance < total_distance)
+      continue;
+    double mfp = p.traversal_mfp();
+    double attenuation = std::exp(-mfp);
+
+    // Save the attenuation for point filter handling
+    p.wgt_last() = p.wgt();
+    p.wgt() *= attenuation;
+
+    double flux = p.wgt_last() * pdf;
+    score_tracklength_tally_general(p, flux, model::active_point_tallies);
+  }
+}
+
+} // namespace openmc
+
+#endif // OPENMC_TALLIES_NEXT_EVENT_SCORING_H
diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h
index 66e6ff764f3..de8495b5196 100644
--- a/include/openmc/tallies/tally.h
+++ b/include/openmc/tallies/tally.h
@@ -3,6 +3,7 @@
 
 #include "openmc/constants.h"
 #include "openmc/memory.h" // for unique_ptr
+#include "openmc/position.h"
 #include "openmc/span.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/trigger.h"
@@ -11,6 +12,7 @@
 #include "openmc/tensor.h"
 #include "pugixml.hpp"
 
+#include <set>
 #include <string>
 #include <unordered_map>
 
@@ -209,11 +211,13 @@ extern vector<int> active_analog_tallies;
 extern vector<int> active_tracklength_tallies;
 extern vector<int> active_timed_tracklength_tallies;
 extern vector<int> active_collision_tallies;
+extern vector<int> active_point_tallies;
 extern vector<int> active_meshsurf_tallies;
 extern vector<int> active_surface_tallies;
 extern vector<int> active_pulse_height_tallies;
 extern vector<int32_t> pulse_height_cells;
 extern vector<double> time_grid;
+extern std::set<Position> active_point_detectors;
 
 } // namespace model
 
diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h
index d1aed28318c..ec806d32181 100644
--- a/include/openmc/tallies/tally_scoring.h
+++ b/include/openmc/tallies/tally_scoring.h
@@ -1,9 +1,12 @@
 #ifndef OPENMC_TALLIES_TALLY_SCORING_H
 #define OPENMC_TALLIES_TALLY_SCORING_H
 
+#include "openmc/nuclide.h"
 #include "openmc/particle.h"
+#include "openmc/ray.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
+#include "openmc/thermal.h"
 
 namespace openmc {
 
@@ -81,6 +84,9 @@ void score_analog_tally_ce(Particle& p);
 //! \param p The particle being tracked
 void score_analog_tally_mg(Particle& p);
 
+void score_tracklength_tally_general(
+  Particle& p, double flux, const vector<int>& tallies);
+
 //! Score tallies using a tracklength estimate of the flux.
 //
 //! This is triggered at every event (surface crossing, lattice crossing, or
diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h
index c06a2ee0dde..a50af6cfd0b 100644
--- a/include/openmc/thermal.h
+++ b/include/openmc/thermal.h
@@ -67,9 +67,11 @@ class ThermalData {
   //! \param[in] mu Scattering cosine with respect to current direction
   //! \param[out] E_out Outgoing energy in [eV]
   //! \param[inout] seed Pseudorandom seed pointer
-  //! \return Probability density for the scattering cosine
+  //! \param[in] is_com Is scattering cosine is given in center of mass
+  //! coordinates \param[in] awr Weight of nucleus in neutron masses \return
+  //! Probability density for the scattering cosine
   double sample_energy_and_pdf(const NuclideMicroXS& micro_xs, double E_in,
-    double mu, double& E_out, uint64_t* seed) const;
+    double mu, double& E_out, uint64_t* seed, bool is_com, double awr) const;
 
 private:
   struct Reaction {
diff --git a/openmc/filter.py b/openmc/filter.py
index 87aeb70c3a7..23f1d6dd16e 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -26,7 +26,7 @@
     'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell',
     'delayedgroup', 'energyfunction', 'cellfrom', 'materialfrom', 'legendre',
     'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle',
-    'particleproduction', 'cellinstance', 'collision', 'time', 'parentnuclide',
+    'particleproduction', 'point', 'cellinstance', 'collision', 'time', 'parentnuclide',
     'weight', 'meshborn', 'meshsurface', 'meshmaterial', 'reaction',
 )
 
@@ -793,6 +793,88 @@ def from_xml_element(cls, elem, **kwargs):
         return cls(bins, filter_id=filter_id)
 
 
+class PointFilter(Filter):
+    """Bins tally events based on point detectors.
+
+    Parameters
+    ----------
+    bins : sequence of tuple[tuple[Real, Real, Real], Real]
+        Point detectors positions and exclusion radii.
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    bins : sequence of tuple[tuple[Real, Real, Real], Real]
+        Point detectors positions and exclusion radii.
+    id : int
+        Unique identifier for the filter
+    num_bins : Integral
+        The number of filter bins
+
+    """
+    
+    __hash__ = Filter.__hash__
+    
+    def __eq__(self, other):
+        if type(self) is not type(other):
+            return False
+        elif len(self.bins) != len(other.bins):
+            return False
+        else:
+            return all(b1==b2 for b1,b2 in zip(self.bins,other.bins))
+
+    @Filter.bins.setter
+    def bins(self, bins):
+        cv.check_type('bins', bins, Sequence, tuple)
+        for i, item in enumerate(bins):
+            cv.check_type(f'bins[{i}]', item, tuple)
+            cv.check_length(f'bins[{i}]', item, 2, 2)
+            cv.check_type(f'bins[{i}][0]', item[0], tuple, Real)
+            cv.check_length(f'bins[{i}][0]', item[0], 3, 3)
+            cv.check_type(f'bins[{i}][1]', item[1], Real)
+        self._bins = bins
+
+    @classmethod
+    def from_hdf5(cls, group, **kwargs):
+        filter_id = int(group.name.split('/')[-1].lstrip('filter '))
+        flat = group['bins'][()]
+        # Reconstruct tuple structure: every 4 values = (x, y, z, r0)
+        bins = []
+        for i in range(0, len(flat), 4):
+            pos = (float(flat[i]), float(flat[i+1]), float(flat[i+2]))
+            r0 = float(flat[i+3])
+            bins.append((pos, r0))
+        out = cls(bins, filter_id=filter_id)
+        out._num_bins = group['n_bins'][()]
+        return out
+
+    def get_pandas_dataframe(self, data_size, stride, **kwargs):
+        import pandas as pd
+        labels = [f"({p[0]}, {p[1]}, {p[2]}) R0={r}" for (p, r) in self.bins]
+        filter_bins = np.repeat(labels, stride)
+        tile_factor = data_size // len(filter_bins)
+        filter_bins = np.tile(filter_bins, tile_factor)
+        return pd.DataFrame({self.short_name.lower(): filter_bins})
+    
+    def to_xml_element(self):
+        """Return XML Element representing the Filter.
+
+        Returns
+        -------
+        element : lxml.etree._Element
+            XML element containing filter data
+
+        """
+        element = ET.Element('filter')
+        element.set('id', str(self.id))
+        element.set('type', self.short_name.lower())
+
+        subelement = ET.SubElement(element, 'bins')
+        subelement.text = ' '.join(str(b) for item in self.bins 
+                                          for b in list(item[0])+[item[1]])
+        return element    
+
 class ParentNuclideFilter(ParticleFilter):
     """Bins tally events based on the parent nuclide
 
diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py
index 574a37443ae..3db32cd59a5 100644
--- a/openmc/lib/filter.py
+++ b/openmc/lib/filter.py
@@ -22,7 +22,7 @@
     'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter',
     'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter',
     'MeshMaterialFilter', 'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter',
-    'ParentNuclideFilter', 'ParticleFilter', 'ParticleProductionFilter', 'PolarFilter',
+    'ParentNuclideFilter', 'ParticleFilter', 'ParticleProductionFilter', 'PointFilter', 'PolarFilter',
     'ReactionFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter',
     'SurfaceFilter', 'TimeFilter', 'UniverseFilter', 'WeightFilter', 'ZernikeFilter',
     'ZernikeRadialFilter', 'filters'
@@ -607,6 +607,8 @@ def bins(self):
             self._index, particle_i.ctypes.data_as(POINTER(c_int32)))
         return [ParticleType(i) for i in particle_i]
 
+class PointFilter(Filter):
+    filter_type = 'point'
 
 class ParticleProductionFilter(Filter):
     filter_type = 'particleproduction'
diff --git a/openmc/tallies.py b/openmc/tallies.py
index ceced425533..07b09edb903 100644
--- a/openmc/tallies.py
+++ b/openmc/tallies.py
@@ -53,7 +53,7 @@
 _FILTER_CLASSES = (Filter, CrossFilter, AggregateFilter)
 
 # Valid types of estimators
-ESTIMATOR_TYPES = {'tracklength', 'collision', 'analog'}
+ESTIMATOR_TYPES = {'tracklength', 'collision', 'analog', 'next-event'}
 
 
 class Tally(IDManagerMixin):
diff --git a/src/angle_energy.cpp b/src/angle_energy.cpp
new file mode 100644
index 00000000000..a3c9c83af77
--- /dev/null
+++ b/src/angle_energy.cpp
@@ -0,0 +1,50 @@
+#include "openmc/angle_energy.h"
+
+#include <algorithm> // for clamp
+#include <cmath>     // for sqrt
+
+#include "openmc/constants.h"
+#include "openmc/random_lcg.h"
+
+namespace openmc {
+
+double get_jac_and_transform(
+  double E_in, double& mu, double& E_out, uint64_t* seed, double awr)
+{
+  double E_com = E_in / ((awr + 1.0) * (awr + 1.0));
+  return get_jac_and_transform_impl(E_com, mu, E_out, seed, awr);
+}
+
+double get_jac_and_transform_impl(
+  double E_com, double& mu, double& E_out, uint64_t* seed, double awr)
+{
+  double E_cm = E_out;
+  double mu_lab = mu;
+  double D = mu_lab * mu_lab - 1.0 + E_cm / E_com;
+  if (D <= 0.0)
+    return 0.0;
+  D = std::sqrt(D);
+
+  if ((mu_lab <= 0.0) && (E_cm <= E_com))
+    return 0.0;
+  double E_out1 = E_com * (mu_lab + D) * (mu_lab + D);
+  double mult;
+  if (E_cm > E_com) {
+    mult = 1.0;
+    E_out = E_out1;
+  } else {
+    mult = 2.0;
+    if (prn(seed) < 0.5) {
+      E_out = E_out1;
+    } else {
+      E_out = E_com * (mu_lab - D) * (mu_lab - D);
+    }
+  }
+  mu = mu_lab * std::sqrt(E_out / E_cm) - std::sqrt(E_com / E_cm);
+
+  if ((std::abs(mu) > 1.0) && (std::abs(mu) < 1.0 + FP_PRECISION))
+    mu = std::clamp(mu, -1.0, 1.0);
+  return mult * E_out / (D * std::sqrt(E_cm * E_com));
+}
+
+} // namespace openmc
diff --git a/src/chain.cpp b/src/chain.cpp
index 3915c016c29..026549208f0 100644
--- a/src/chain.cpp
+++ b/src/chain.cpp
@@ -74,9 +74,13 @@ void DecayPhotonAngleEnergy::sample(
   mu = Uniform(-1., 1.).sample(seed).first;
 }
 
-double DecayPhotonAngleEnergy::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double DecayPhotonAngleEnergy::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
+  if (is_com)
+    fatal_error(
+      "DecayPhotonAngleEnergy distribution must be in lab coordinates");
+
   E_out = photon_energy_->sample(seed).first;
   return 0.5;
 }
diff --git a/src/physics.cpp b/src/physics.cpp
index 106bd1aa2b2..fb686c540a3 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -24,7 +24,9 @@
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
 #include "openmc/string_utils.h"
+#include "openmc/tallies/next_event_scoring.h"
 #include "openmc/tallies/tally.h"
+#include "openmc/tallies/tally_scoring.h"
 #include "openmc/thermal.h"
 #include "openmc/weight_windows.h"
 
@@ -739,7 +741,7 @@ void scatter(Particle& p, int i_nuclide)
 
     // Perform collision physics for inelastic scattering
     const auto& rx {nuc->reactions_[i]};
-    inelastic_scatter(*nuc, *rx, p);
+    inelastic_scatter(i_nuclide, *rx, p);
     p.event_mt() = rx->mt_;
   }
 
@@ -785,6 +787,10 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, Particle& p)
   // Find speed of neutron in CM
   vel = v_n.norm();
 
+  if (!model::active_point_tallies.empty()) {
+    score_point_tally_elastic(p, i_nuclide, rx, 0, v_t);
+  }
+
   // Sample scattering angle, checking if angle distribution is present (assume
   // isotropic otherwise)
   double mu_cm;
@@ -832,8 +838,13 @@ void sab_scatter(int i_nuclide, int i_sab, Particle& p)
 
   // Sample energy and angle
   double E_out;
-  data::thermal_scatt[i_sab]->data_[i_temp].sample(
-    micro, p.E(), &E_out, &p.mu(), p.current_seed());
+  auto& sab = data::thermal_scatt[i_sab]->data_[i_temp];
+
+  if (!model::active_point_tallies.empty()) {
+    score_point_tally_sab(p, i_nuclide, sab, micro);
+  }
+
+  sab.sample(micro, p.E(), &E_out, &p.mu(), p.current_seed());
 
   // Set energy to outgoing, change direction of particle
   p.E() = E_out;
@@ -1098,6 +1109,10 @@ void sample_fission_neutron(
     site->delayed_group = 0;
   }
 
+  if (!model::active_point_tallies.empty()) {
+    score_point_tally_fission(p, i_nuclide, rx, site->delayed_group);
+  }
+
   // sample from prompt neutron energy distribution
   int n_sample = 0;
   double mu;
@@ -1123,8 +1138,11 @@ void sample_fission_neutron(
   site->u = rotate_angle(p.u(), mu, nullptr, seed);
 }
 
-void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
+void inelastic_scatter(int i_nuclide, const Reaction& rx, Particle& p)
 {
+  // Get pointer to nuclide
+  const auto& nuc {data::nuclides[i_nuclide]};
+
   // copy energy of neutron
   double E_in = p.E();
 
@@ -1133,13 +1151,19 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
   double mu;
   rx.products_[0].sample(E_in, E, mu, p.current_seed());
 
+  double yield = (*rx.products_[0].yield_)(E_in);
+
+  if (!model::active_point_tallies.empty()) {
+    score_point_tally_inelastic(p, i_nuclide, rx, 0, yield);
+  }
+
   // if scattering system is in center-of-mass, transfer cosine of scattering
   // angle and outgoing energy from CM to LAB
   if (rx.scatter_in_cm_) {
     double E_cm = E;
 
     // determine outgoing energy in lab
-    double A = nuc.awr_;
+    double A = nuc->awr_;
     E = E_cm + (E_in + 2.0 * mu * (A + 1.0) * std::sqrt(E_in * E_cm)) /
                  ((A + 1.0) * (A + 1.0));
 
@@ -1160,8 +1184,6 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
   // change direction of particle
   p.u() = rotate_angle(p.u(), mu, nullptr, p.current_seed());
 
-  // evaluate yield
-  double yield = (*rx.products_[0].yield_)(E_in);
   if (std::floor(yield) == yield && yield > 0) {
     // If yield is integral, create exactly that many secondary particles
     for (int i = 0; i < static_cast<int>(std::round(yield)) - 1; ++i) {
diff --git a/src/ray.cpp b/src/ray.cpp
index 3d848e3a3a7..5e9a8419ff7 100644
--- a/src/ray.cpp
+++ b/src/ray.cpp
@@ -2,6 +2,8 @@
 
 #include "openmc/error.h"
 #include "openmc/geometry.h"
+#include "openmc/material.h"
+#include "openmc/mgxs_interface.h"
 #include "openmc/settings.h"
 
 namespace openmc {
@@ -11,7 +13,7 @@ void Ray::compute_distance()
   boundary() = distance_to_boundary(*this);
 }
 
-void Ray::trace()
+void Ray::trace(double max_distance)
 {
   // To trace the ray from its origin all the way through the model, we have
   // to proceed in two phases. In the first, the ray may or may not be found
@@ -24,6 +26,8 @@ void Ray::trace()
   // After phase one is done, we can starting tracing from cell to cell within
   // the model. This step can use neighbor lists to accelerate the ray tracing.
 
+  double max = max_distance;
+
   bool inside_cell;
   // Check for location if the particle is already known
   if (lowest_coord().cell() == C_NONE) {
@@ -118,10 +122,20 @@ void Ray::trace()
     // distance to properly check cell inclusion.
     boundary().distance() += TINY_BIT;
 
+    double distance = std::min(boundary().distance(), max);
+
     // Advance particle, prepare for next intersection
     for (int lev = 0; lev < n_coord(); ++lev) {
-      coord(lev).r() += boundary().distance() * coord(lev).u();
+      coord(lev).r() += distance * coord(lev).u();
+    }
+
+    max -= distance;
+
+    if (max == 0.0) {
+      update_distance();
+      break;
     }
+
     surface() = boundary().surface();
     // Initialize last cells from the current cell, because the cell() variable
     // does not contain the data for the case of a single-segment ray
@@ -136,8 +150,8 @@ void Ray::trace()
       cross_lattice(*this, boundary(), settings::verbosity >= 10);
     }
 
-    // Record how far the ray has traveled
-    traversal_distance_ += boundary().distance();
+    update_distance();
+
     inside_cell = neighbor_list_find_cell(*this, settings::verbosity >= 10);
 
     // Call the specialized logic for this type of ray. Note that we do not
@@ -165,4 +179,47 @@ void Ray::trace()
   }
 }
 
+void Ray::update_distance()
+{
+  // Record how far the ray has traveled
+  traversal_distance_ += boundary().distance();
+}
+
+void ParticleRay::on_intersection() {}
+
+void ParticleRay::update_distance()
+{
+  Ray::update_distance();
+
+  time() += boundary().distance() / speed();
+
+  // Calculate microscopic and macroscopic cross sections
+  if (material() != MATERIAL_VOID) {
+    if (settings::run_CE) {
+      if (material() != material_last() || sqrtkT() != sqrtkT_last() ||
+          density_mult() != density_mult_last()) {
+        // If the material is the same as the last material and the
+        // temperature hasn't changed, we don't need to lookup cross
+        // sections again.
+        model::materials[material()]->calculate_xs(*this);
+      }
+    } else {
+      // Get the MG data; unlike the CE case above, we have to re-calculate
+      // cross sections for every collision since the cross sections may
+      // be angle-dependent
+      data::mg.macro_xs_[material()].calculate_xs(*this);
+
+      // Update the particle's group while we know we are multi-group
+      g_last() = g();
+    }
+  } else {
+    macro_xs().total = 0.0;
+    macro_xs().absorption = 0.0;
+    macro_xs().fission = 0.0;
+    macro_xs().nu_fission = 0.0;
+  }
+
+  traversal_mfp_ += macro_xs().total * boundary().distance();
+}
+
 } // namespace openmc
diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp
index a1c93786166..69fc959b944 100644
--- a/src/reaction_product.cpp
+++ b/src/reaction_product.cpp
@@ -134,10 +134,11 @@ void ReactionProduct::sample(
   sample_dist(E_in, seed).sample(E_in, E_out, mu, seed);
 }
 
-double ReactionProduct::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double ReactionProduct::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
-  return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed);
+  return sample_dist(E_in, seed)
+    .sample_energy_and_pdf(E_in, mu, E_out, seed, is_com, awr);
 }
 
 } // namespace openmc
diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp
index 32701791ade..99b49cb1970 100644
--- a/src/secondary_correlated.cpp
+++ b/src/secondary_correlated.cpp
@@ -260,10 +260,15 @@ void CorrelatedAngleEnergy::sample(
   mu = sample_dist(E_in, E_out, seed).sample(seed).first;
 }
 
-double CorrelatedAngleEnergy::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double CorrelatedAngleEnergy::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
-  return sample_dist(E_in, E_out, seed).evaluate(mu);
+  auto& dist = sample_dist(E_in, E_out, seed);
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
+  return jac * dist.evaluate(mu);
 }
 
 } // namespace openmc
diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp
index 018ce1c8a9c..3472d9e7391 100644
--- a/src/secondary_kalbach.cpp
+++ b/src/secondary_kalbach.cpp
@@ -232,14 +232,18 @@ void KalbachMann::sample(
     mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a;
   }
 }
-double KalbachMann::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double KalbachMann::sample_energy_and_pdf(double E_in, double mu, double& E_out,
+  uint64_t* seed, bool is_com, double awr) const
 {
   double km_r, km_a;
   sample_params(E_in, E_out, km_a, km_r, seed);
 
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
   // https://docs.openmc.org/en/latest/methods/neutron_physics.html#equation-KM-pdf-angle
-  return km_a / (2 * std::sinh(km_a)) *
+  return jac * km_a / (2 * std::sinh(km_a)) *
          (std::cosh(km_a * mu) + km_r * std::sinh(km_a * mu));
 }
 
diff --git a/src/secondary_nbody.cpp b/src/secondary_nbody.cpp
index 72f0b0b92d0..9ae7a5152f4 100644
--- a/src/secondary_nbody.cpp
+++ b/src/secondary_nbody.cpp
@@ -72,11 +72,14 @@ void NBodyPhaseSpace::sample(
   E_out = sample_energy(E_in, seed);
 }
 
-double NBodyPhaseSpace::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double NBodyPhaseSpace::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
   E_out = sample_energy(E_in, seed);
-  return 0.5;
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+  return jac * 0.5;
 }
 
 } // namespace openmc
diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp
index b0f601809df..7634c15ed16 100644
--- a/src/secondary_thermal.cpp
+++ b/src/secondary_thermal.cpp
@@ -53,11 +53,16 @@ void CoherentElasticAE::sample(
   mu = 1.0 - 2.0 * energies[k] / E_in;
 }
 
-double CoherentElasticAE::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double CoherentElasticAE::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
   // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
   E_out = E_in;
+
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
   const auto& factors = xs_.factors();
 
   if (E_in < bragg_edges_.front())
@@ -68,8 +73,9 @@ double CoherentElasticAE::sample_energy_and_pdf(
   double E = 0.5 * (1 - mu) * E_in;
   double C = 0.5 * E_in / factors[i];
 
-  return C * get_pdf_discrete(bragg_edges_.slice(tensor::range(i + 1)),
-               factors_diff_.slice(tensor::range(i + 1)), E, 0.0, E_in);
+  return jac * C *
+         get_pdf_discrete(bragg_edges_.slice(tensor::range(i + 1)),
+           factors_diff_.slice(tensor::range(i + 1)), E, 0.0, E_in);
 }
 
 //==============================================================================
@@ -90,15 +96,20 @@ void IncoherentElasticAE::sample(
   double c = 2 * E_in * debye_waller_;
   mu = std::log(1.0 + prn(seed) * (std::exp(2.0 * c) - 1)) / c - 1.0;
 }
-double IncoherentElasticAE::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+
+double IncoherentElasticAE::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
   E_out = E_in;
 
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
   // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
   double c = 2 * E_in * debye_waller_;
   double A = c / (1 - std::exp(-2.0 * c)); // normalization factor
-  return A * std::exp(-c * (1 - mu));
+  return jac * A * std::exp(-c * (1 - mu));
 }
 
 //==============================================================================
@@ -155,8 +166,8 @@ void IncoherentElasticAEDiscrete::sample(
   E_out = E_in;
 }
 
-double IncoherentElasticAEDiscrete::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double IncoherentElasticAEDiscrete::sample_energy_and_pdf(double E_in,
+  double mu, double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
   // Get index and interpolation factor for elastic grid
   int i;
@@ -165,8 +176,12 @@ double IncoherentElasticAEDiscrete::sample_energy_and_pdf(
   // Energy doesn't change in elastic scattering
   E_out = E_in;
 
-  return get_pdf_discrete_interpolated(
-    mu_out_.slice(i, tensor::all), mu_out_.slice(i + 1, tensor::all), f, mu);
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
+  return jac * get_pdf_discrete_interpolated(mu_out_.slice(i, tensor::all),
+                 mu_out_.slice(i + 1, tensor::all), f, mu);
 }
 
 //==============================================================================
@@ -253,8 +268,8 @@ void IncoherentInelasticAEDiscrete::sample(
   mu = (1 - f) * mu_ijk + f * mu_i1jk;
 }
 
-double IncoherentInelasticAEDiscrete::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double IncoherentInelasticAEDiscrete::sample_energy_and_pdf(double E_in,
+  double mu, double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
   // Get index and interpolation factor for inelastic grid
   int i;
@@ -263,8 +278,12 @@ double IncoherentInelasticAEDiscrete::sample_energy_and_pdf(
   int j;
   sample_params(E_in, E_out, j, seed);
 
-  return get_pdf_discrete_interpolated(mu_out_.slice(i, j, tensor::all),
-    mu_out_.slice(i + 1, j, tensor::all), f, mu);
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
+  return jac * get_pdf_discrete_interpolated(mu_out_.slice(i, j, tensor::all),
+                 mu_out_.slice(i + 1, j, tensor::all), f, mu);
 }
 
 //==============================================================================
@@ -403,8 +422,8 @@ void IncoherentInelasticAE::sample(
   mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5);
 }
 
-double IncoherentInelasticAE::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double IncoherentInelasticAE::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
   double f;
   int l, j;
@@ -412,8 +431,12 @@ double IncoherentInelasticAE::sample_energy_and_pdf(
 
   const auto& mu_l = distribution_[l].mu;
 
-  return get_pdf_discrete_interpolated(
-    mu_l.slice(j, tensor::all), mu_l.slice(j + 1, tensor::all), f, mu);
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
+  return jac * get_pdf_discrete_interpolated(mu_l.slice(j, tensor::all),
+                 mu_l.slice(j + 1, tensor::all), f, mu);
 }
 
 //==============================================================================
@@ -458,10 +481,11 @@ void MixedElasticAE::sample(
   sample_dist(E_in, seed).sample(E_in, E_out, mu, seed);
 }
 
-double MixedElasticAE::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double MixedElasticAE::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
-  return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed);
+  return sample_dist(E_in, seed)
+    .sample_energy_and_pdf(E_in, mu, E_out, seed, is_com, awr);
 }
 
 } // namespace openmc
diff --git a/src/secondary_uncorrelated.cpp b/src/secondary_uncorrelated.cpp
index ec2af710258..59818832247 100644
--- a/src/secondary_uncorrelated.cpp
+++ b/src/secondary_uncorrelated.cpp
@@ -65,8 +65,8 @@ void UncorrelatedAngleEnergy::sample(
   E_out = energy_->sample(E_in, seed);
 }
 
-double UncorrelatedAngleEnergy::sample_energy_and_pdf(
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+double UncorrelatedAngleEnergy::sample_energy_and_pdf(double E_in, double mu,
+  double& E_out, uint64_t* seed, bool is_com, double awr) const
 {
   // Sample outgoing energy
   if (energy_ != nullptr) {
@@ -75,11 +75,15 @@ double UncorrelatedAngleEnergy::sample_energy_and_pdf(
     E_out = E_in;
   }
 
+  double jac = 1.0;
+  if (is_com)
+    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
+
   if (!angle_.empty()) {
-    return angle_.evaluate(E_in, mu);
+    return jac * angle_.evaluate(E_in, mu);
   } else {
     // no angle distribution given => assume isotropic for all energies
-    return 0.5;
+    return jac * 0.5;
   }
 }
 
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 4fad196a604..c384b86f04c 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -311,6 +311,7 @@ double k_abs_tra {0.0};
 double log_spacing;
 int n_lost_particles {0};
 bool need_depletion_rx {false};
+bool nonvacuum_boundary_present {false};
 int restart_batch;
 bool satisfy_triggers {false};
 int ssw_current_file;
diff --git a/src/source.cpp b/src/source.cpp
index bee20d5c305..dba065cd5d7 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -33,6 +33,8 @@
 #include "openmc/simulation.h"
 #include "openmc/state_point.h"
 #include "openmc/string_utils.h"
+#include "openmc/tallies/next_event_scoring.h"
+#include "openmc/tallies/tally_scoring.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -728,6 +730,10 @@ SourceSite sample_external_source(uint64_t* seed)
     site.E = data::mg.num_energy_groups_ - site.E - 1.;
   }
 
+  if (!model::active_point_tallies.empty()) {
+    score_point_tally_source(site, i);
+  }
+
   return site;
 }
 
diff --git a/src/state_point.cpp b/src/state_point.cpp
index da1c141a230..f599ada7947 100644
--- a/src/state_point.cpp
+++ b/src/state_point.cpp
@@ -214,6 +214,8 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
           write_dataset(tally_group, "estimator", "tracklength");
         } else if (tally->estimator_ == TallyEstimator::COLLISION) {
           write_dataset(tally_group, "estimator", "collision");
+        } else if (tally->estimator_ == TallyEstimator::NEXT_EVENT) {
+          write_dataset(tally_group, "estimator", "next-event");
         }
 
         write_dataset(tally_group, "n_realizations", tally->n_realizations_);
diff --git a/src/surface.cpp b/src/surface.cpp
index 81b756deae7..08cdc494807 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -17,6 +17,7 @@
 #include "openmc/math_functions.h"
 #include "openmc/random_lcg.h"
 #include "openmc/settings.h"
+#include "openmc/simulation.h"
 #include "openmc/string_utils.h"
 #include "openmc/xml_interface.h"
 
@@ -1175,6 +1176,7 @@ void read_surfaces(pugi::xml_node node,
   std::unordered_map<int, double>& albedo_map,
   std::unordered_map<int, int>& periodic_sense_map)
 {
+  simulation::nonvacuum_boundary_present = false;
   // Count the number of surfaces
   int n_surfaces = 0;
   for (pugi::xml_node surf_node : node.children("surface")) {
@@ -1245,6 +1247,8 @@ void read_surfaces(pugi::xml_node node,
       // Check for a periodic surface
       if (check_for_node(surf_node, "boundary")) {
         std::string surf_bc = get_node_value(surf_node, "boundary", true, true);
+        if (surf_bc != "vacuum")
+          simulation::nonvacuum_boundary_present = true;
         if (surf_bc == "periodic") {
           periodic_sense_map[model::surfaces.back()->id_] = 0;
           // Check for surface albedo. Skip sanity check as it is already done
diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp
index badb9107733..1068b4b5d08 100644
--- a/src/tallies/filter.cpp
+++ b/src/tallies/filter.cpp
@@ -32,6 +32,7 @@
 #include "openmc/tallies/filter_parent_nuclide.h"
 #include "openmc/tallies/filter_particle.h"
 #include "openmc/tallies/filter_particle_production.h"
+#include "openmc/tallies/filter_point.h"
 #include "openmc/tallies/filter_polar.h"
 #include "openmc/tallies/filter_reaction.h"
 #include "openmc/tallies/filter_sph_harm.h"
@@ -150,6 +151,8 @@ Filter* Filter::create(const std::string& type, int32_t id)
     return Filter::create<ParticleFilter>(id);
   } else if (type == "particleproduction") {
     return Filter::create<ParticleProductionFilter>(id);
+  } else if (type == "point") {
+    return Filter::create<PointFilter>(id);
   } else if (type == "polar") {
     return Filter::create<PolarFilter>(id);
   } else if (type == "reaction") {
diff --git a/src/tallies/filter_point.cpp b/src/tallies/filter_point.cpp
new file mode 100644
index 00000000000..35f4d2e03cf
--- /dev/null
+++ b/src/tallies/filter_point.cpp
@@ -0,0 +1,79 @@
+#include "openmc/tallies/filter_point.h"
+#include "openmc/tallies/tally_scoring.h"
+
+#include <fmt/core.h>
+
+#include "openmc/math_functions.h"
+#include "openmc/simulation.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void PointFilter::from_xml(pugi::xml_node node)
+{
+  auto bins = get_node_array<double>(node, "bins");
+
+  // Convert to vector of detectors
+  vector<std::pair<Position, double>> detectors;
+  size_t n = bins.size() / 4;
+  for (int i = 0; i < n; ++i) {
+    Position pos {bins[4 * i], bins[4 * i + 1], bins[4 * i + 2]};
+    detectors.push_back(std::make_pair(pos, bins[4 * i + 3]));
+  }
+  this->set_detectors(detectors);
+}
+
+void PointFilter::set_detectors(span<std::pair<Position, double>> detectors)
+{
+  // Clear existing detectors
+  detectors_.clear();
+  detectors_.reserve(detectors.size());
+
+  // Set detectors and number of bins
+  for (auto d : detectors) {
+    detectors_.push_back(d);
+  }
+  n_bins_ = detectors_.size();
+}
+
+void PointFilter::get_all_bins(
+  const Particle& p, TallyEstimator estimator, FilterMatch& match) const
+{
+  double attenuation = p.wgt() / p.wgt_last();
+  int i = 0;
+  for (auto [pos, r] : detectors_) {
+    if ((p.r() - pos).norm() < FP_COINCIDENT) {
+      match.bins_.push_back(i);
+      double weight;
+      double distance = (p.r_last() - pos).norm();
+      if (distance > r) {
+        weight = attenuation / (distance * distance);
+      } else {
+        weight = 3.0 * exprel(-p.macro_xs().total * r) / (r * r);
+      }
+      match.weights_.push_back(weight);
+    }
+    ++i;
+  }
+}
+
+void PointFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  vector<double> detectors;
+  for (auto [pos, r] : detectors_) {
+    detectors.push_back(pos[0]);
+    detectors.push_back(pos[1]);
+    detectors.push_back(pos[2]);
+    detectors.push_back(r);
+  }
+  write_dataset(filter_group, "bins", detectors);
+}
+
+std::string PointFilter::text_label(int bin) const
+{
+  auto [pos, r] = detectors_.at(bin);
+  return fmt::format("Point: {} {}", pos, r);
+}
+
+} // namespace openmc
diff --git a/src/tallies/next_event_scoring.cpp b/src/tallies/next_event_scoring.cpp
new file mode 100644
index 00000000000..26ea2c18019
--- /dev/null
+++ b/src/tallies/next_event_scoring.cpp
@@ -0,0 +1,115 @@
+#include "openmc/tallies/next_event_scoring.h"
+
+#include <algorithm> // for clamp
+
+#include "openmc/secondary_uncorrelated.h"
+#include "openmc/source.h"
+
+namespace openmc {
+
+void score_point_tally_elastic(
+  Particle& p, int i_nuclide, const Reaction& rx, int i_product, Direction v_t)
+{
+
+  const auto& nuc {data::nuclides[i_nuclide]};
+  double awr = nuc->awr_;
+
+  // Neutron velocity in LAB
+  Direction v_n = std::sqrt(p.E()) * p.u();
+  auto u_n = v_n / v_n.norm();
+
+  // Velocity of center-of-mass
+  Direction v_cm = (v_n + awr * v_t) / (awr + 1.0);
+  auto u_cm = v_cm / v_cm.norm();
+
+  double E_in = p.E();
+  double E_com = v_cm.dot(v_cm);
+  double E_out = (v_n - v_cm).dot(v_n - v_cm);
+
+  auto& d = rx.products_[i_product].distribution_[0];
+  auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
+
+  auto pdf = [&](Direction u, double& E) {
+    double mu = u.dot(u_cm);
+    double mu_l = u.dot(u_n);
+    E = E_out;
+    double jac =
+      get_jac_and_transform_impl(E_com, mu, E, p.current_seed(), awr);
+    double mu_cm =
+      1.0 + mu_l * std::sqrt(E_in * E) / E_out - (E_in + E) / (2.0 * E_out);
+    mu_cm = std::clamp(mu_cm, -1.0, 1.0);
+    if (!d_->angle().empty()) {
+      return jac * d_->angle().evaluate(p.E(), mu_cm) / (2.0 * PI);
+    } else {
+      return jac * 0.5 / (2.0 * PI);
+    }
+  };
+  score_point_tally_impl(p.r(), p.type(), p.time(), pdf);
+}
+
+void score_point_tally_inelastic(
+  Particle& p, int i_nuclide, const Reaction& rx, int i_product, double yield)
+{
+  const auto& nuc {data::nuclides[i_nuclide]};
+  double awr = nuc->awr_;
+  auto u_n = p.u();
+  auto E_n = p.E();
+  auto is_com = rx.scatter_in_cm_;
+
+  auto pdf = [&](Direction u, double& E) {
+    double mu = u.dot(u_n);
+    return rx.products_[i_product].sample_energy_and_pdf(
+             E_n, mu, E, p.current_seed(), is_com, awr) /
+           (2.0 * PI) * yield;
+  };
+  score_point_tally_impl(p.r(), p.type(), p.time(), pdf);
+}
+
+void score_point_tally_fission(
+  Particle& p, int i_nuclide, const Reaction& rx, int i_product)
+{
+  const auto& nuc {data::nuclides[i_nuclide]};
+  double awr = nuc->awr_;
+  auto u_n = p.u();
+  auto E_n = p.E();
+  auto is_com = rx.scatter_in_cm_;
+
+  auto pdf = [&](Direction u, double& E) {
+    double mu = u.dot(u_n);
+    return rx.products_[i_product].sample_energy_and_pdf(
+             E_n, mu, E, p.current_seed(), is_com, awr) /
+           (2.0 * PI);
+  };
+  score_point_tally_impl(p.r(), p.type(), p.time(), pdf);
+}
+
+void score_point_tally_sab(Particle& p, int i_nuclide, const ThermalData& sab,
+  const NuclideMicroXS& micro)
+{
+  const auto& nuc {data::nuclides[i_nuclide]};
+  double awr = nuc->awr_;
+  auto u_n = p.u();
+  auto E_n = p.E();
+  auto pdf = [&](Direction u, double& E) {
+    double mu = u.dot(u_n);
+    return sab.sample_energy_and_pdf(
+             micro, E_n, mu, E, p.current_seed(), false, awr) /
+           (2.0 * PI);
+  };
+  score_point_tally_impl(p.r(), p.type(), p.time(), pdf);
+}
+
+void score_point_tally_source(SourceSite& site, int source_index)
+{
+  auto src_ = model::external_sources[source_index].get();
+  auto src = dynamic_cast<IndependentSource*>(src_);
+  if (!src)
+    fatal_error("Only independent source is valid for point detectors.");
+  auto pdf = [&](Direction u, double& E) {
+    E = site.E;
+    return src->angle()->evaluate(u);
+  };
+  score_point_tally_impl(site.r, site.particle, site.time, pdf);
+}
+
+} // namespace openmc
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index 3fe48c1b021..f4de1253500 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -31,6 +31,7 @@
 #include "openmc/tallies/filter_meshmaterial.h"
 #include "openmc/tallies/filter_meshsurface.h"
 #include "openmc/tallies/filter_particle.h"
+#include "openmc/tallies/filter_point.h"
 #include "openmc/tallies/filter_sph_harm.h"
 #include "openmc/tallies/filter_surface.h"
 #include "openmc/tallies/filter_time.h"
@@ -60,11 +61,13 @@ vector<int> active_analog_tallies;
 vector<int> active_tracklength_tallies;
 vector<int> active_timed_tracklength_tallies;
 vector<int> active_collision_tallies;
+vector<int> active_point_tallies;
 vector<int> active_meshsurf_tallies;
 vector<int> active_surface_tallies;
 vector<int> active_pulse_height_tallies;
 vector<int32_t> pulse_height_cells;
 vector<double> time_grid;
+std::set<Position> active_point_detectors;
 } // namespace model
 
 namespace simulation {
@@ -540,6 +543,7 @@ void Tally::set_scores(const vector<std::string>& scores)
   bool legendre_present = false;
   bool cell_present = false;
   bool cellfrom_present = false;
+  bool point_present = false;
   bool material_present = false;
   bool materialfrom_present = false;
   bool surface_present = false;
@@ -562,6 +566,10 @@ void Tally::set_scores(const vector<std::string>& scores)
       materialfrom_present = true;
     } else if (filt->type() == FilterType::MATERIAL) {
       material_present = true;
+    } else if (filt->type() == FilterType::POINT) {
+      point_present = true;
+      type_ = TallyType::POINT;
+      estimator_ = TallyEstimator::NEXT_EVENT;
     } else if (filt->type() == FilterType::SURFACE) {
       surface_present = true;
     } else if (filt->type() == FilterType::MESH_SURFACE) {
@@ -574,6 +582,19 @@ void Tally::set_scores(const vector<std::string>& scores)
     (surface_present || cell_present || cellfrom_present || material_present ||
       materialfrom_present);
 
+  if (point_present) {
+    if (simulation::nonvacuum_boundary_present)
+      fatal_error(
+        "Cannot use point detectors with non-vacuum boundary conditions.");
+    if (legendre_present)
+      fatal_error("Cannot use LegendreFilter with PointFilter.");
+    if (energyout_present)
+      fatal_error("Cannot use EnergyoutFilter with PointFilter.");
+    if (surface_present || meshsurface_present)
+      fatal_error(
+        "Cannot use surface or mesh-surface filters with PointFilter.");
+  }
+
   // Iterate over the given scores.
   for (auto score_str : scores) {
     // Make sure a delayed group filter wasn't used with an incompatible
@@ -623,6 +644,9 @@ void Tally::set_scores(const vector<std::string>& scores)
 
     case SCORE_NU_SCATTER:
       if (settings::run_CE) {
+        if (point_present)
+          fatal_error("Cannot use nu-scatter score with PointFilter in "
+                      "continuous energy mode.");
         estimator_ = TallyEstimator::ANALOG;
       } else {
         if (energyout_present || legendre_present)
@@ -631,6 +655,8 @@ void Tally::set_scores(const vector<std::string>& scores)
       break;
 
     case SCORE_CURRENT:
+      if (point_present)
+        fatal_error("Cannot use current score with PointFilter.");
       // Check which type of current is desired: mesh or surface currents.
       if (meshsurface_present) {
         if (non_meshsurface_types_present)
@@ -644,11 +670,15 @@ void Tally::set_scores(const vector<std::string>& scores)
       break;
 
     case HEATING:
+      if (point_present)
+        fatal_error("Cannot use heating score with PointFilter.");
       if (settings::photon_transport)
         estimator_ = TallyEstimator::COLLISION;
       break;
 
     case SCORE_PULSE_HEIGHT: {
+      if (point_present)
+        fatal_error("Cannot use pulse-height score with PointFilter.");
       if (non_cell_energy_present) {
         fatal_error("Pulse-height tallies are not compatible with filters "
                     "other than CellFilter and EnergyFilter");
@@ -670,6 +700,8 @@ void Tally::set_scores(const vector<std::string>& scores)
     case SCORE_IFP_TIME_NUM:
     case SCORE_IFP_BETA_NUM:
     case SCORE_IFP_DENOM:
+      if (point_present)
+        fatal_error("Cannot use ifp scores with PointFilter.");
       estimator_ = TallyEstimator::COLLISION;
       break;
     }
@@ -1168,6 +1200,8 @@ void setup_active_tallies()
   model::active_meshsurf_tallies.clear();
   model::active_surface_tallies.clear();
   model::active_pulse_height_tallies.clear();
+  model::active_point_tallies.clear();
+  model::active_point_detectors.clear();
   model::time_grid.clear();
 
   for (auto i = 0; i < model::tallies.size(); ++i) {
@@ -1209,6 +1243,16 @@ void setup_active_tallies()
       case TallyType::PULSE_HEIGHT:
         model::active_pulse_height_tallies.push_back(i);
         break;
+
+      case TallyType::POINT:
+        model::active_point_tallies.push_back(i);
+        // Populate the set of unique detector positions from PointFilter
+        if (auto pf = tally.get_filter<PointFilter>()) {
+          for (const auto& [pos, r0] : pf->detectors()) {
+            model::active_point_detectors.insert(pos);
+          }
+        }
+        break;
       }
     }
   }
@@ -1232,6 +1276,8 @@ void free_memory_tally()
   model::active_meshsurf_tallies.clear();
   model::active_surface_tallies.clear();
   model::active_pulse_height_tallies.clear();
+  model::active_point_tallies.clear();
+  model::active_point_detectors.clear();
   model::time_grid.clear();
 
   model::tally_map.clear();
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index 5dcd07331d6..2c08a7d67dc 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -11,8 +11,10 @@
 #include "openmc/photon.h"
 #include "openmc/reaction_product.h"
 #include "openmc/search.h"
+#include "openmc/secondary_uncorrelated.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
+#include "openmc/source.h"
 #include "openmc/string_utils.h"
 #include "openmc/surface.h"
 #include "openmc/tallies/derivative.h"
@@ -2784,4 +2786,5 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
     p.E_last() = orig_E_last;
   }
 }
+
 } // namespace openmc
diff --git a/src/thermal.cpp b/src/thermal.cpp
index edfbddf23e8..fd753165a39 100644
--- a/src/thermal.cpp
+++ b/src/thermal.cpp
@@ -314,10 +314,11 @@ void ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
 }
 
 double ThermalData::sample_energy_and_pdf(const NuclideMicroXS& micro_xs,
-  double E_in, double mu, double& E_out, uint64_t* seed) const
+  double E_in, double mu, double& E_out, uint64_t* seed, bool is_com,
+  double awr) const
 {
   return sample_dist(micro_xs, E_in, seed)
-    .sample_energy_and_pdf(E_in, mu, E_out, seed);
+    .sample_energy_and_pdf(E_in, mu, E_out, seed, is_com, awr);
 }
 
 void free_memory_thermal()