diff --git a/CMakeLists.txt b/CMakeLists.txt
index 87b8789d101..bd9e93e7d7c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -425,6 +425,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_polar.cpp
   src/tallies/filter_sph_harm.cpp
   src/tallies/filter_sptl_legendre.cpp
diff --git a/include/openmc/angle_energy.h b/include/openmc/angle_energy.h
index ac931b1b533..35e5b68b64b 100644
--- a/include/openmc/angle_energy.h
+++ b/include/openmc/angle_energy.h
@@ -16,6 +16,8 @@ class AngleEnergy {
 public:
   virtual void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const = 0;
+  virtual double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const = 0;
   virtual ~AngleEnergy() = default;
 };
 
diff --git a/include/openmc/chain.h b/include/openmc/chain.h
index a3bc6f3a364..c2a45f75adc 100644
--- a/include/openmc/chain.h
+++ b/include/openmc/chain.h
@@ -70,6 +70,8 @@ class DecayPhotonAngleEnergy : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const Distribution* photon_energy_;
diff --git a/include/openmc/constants.h b/include/openmc/constants.h
index bba260c3a4b..32bf356cb57 100644
--- a/include/openmc/constants.h
+++ b/include/openmc/constants.h
@@ -294,9 +294,15 @@ 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,
+  NEXT_EVENT
+};
 
-enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION };
+enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION, POINT };
 
 enum class TallyEvent { SURFACE, LATTICE, KILL, SCATTER, ABSORB };
 
diff --git a/include/openmc/distribution_angle.h b/include/openmc/distribution_angle.h
index efd4e58425b..a7e021085f0 100644
--- a/include/openmc/distribution_angle.h
+++ b/include/openmc/distribution_angle.h
@@ -25,6 +25,7 @@ class AngleDistribution {
   //! \param[inout] seed pseudorandom number seed pointer
   //! \return Cosine of the angle in the range [-1,1]
   double sample(double E, uint64_t* seed) const;
+  double evaluate(double E, double mu) const;
 
   //! Determine whether angle distribution is empty
   //! \return Whether distribution is empty
diff --git a/include/openmc/distribution_multi.h b/include/openmc/distribution_multi.h
index 7b9c2abf8ce..d53e5ed560e 100644
--- a/include/openmc/distribution_multi.h
+++ b/include/openmc/distribution_multi.h
@@ -6,6 +6,7 @@
 #include "pugixml.hpp"
 
 #include "openmc/distribution.h"
+#include "openmc/error.h"
 #include "openmc/position.h"
 
 namespace openmc {
@@ -29,6 +30,11 @@ class UnitSphereDistribution {
   //! \return (sampled Direction, sample weight)
   virtual std::pair<Direction, double> sample(uint64_t* seed) const = 0;
 
+  virtual double evaluate(Direction u) const
+  {
+    fatal_error("evaluate not available for this UnitSphereDistribution type");
+  }
+
   Direction u_ref_ {0.0, 0.0, 1.0}; //!< reference direction
 };
 
@@ -52,6 +58,8 @@ class PolarAzimuthal : public UnitSphereDistribution {
   //! \return (sampled Direction, value of the PDF at this Direction)
   std::pair<Direction, double> sample_as_bias(uint64_t* seed) const;
 
+  double evaluate(Direction u) const override;
+
   // Observing pointers
   Distribution* mu() const { return mu_.get(); }
   Distribution* phi() const { return phi_.get(); }
@@ -87,6 +95,8 @@ class Isotropic : public UnitSphereDistribution {
   //! \return (sampled direction, sample weight)
   std::pair<Direction, double> sample(uint64_t* seed) const override;
 
+  double evaluate(Direction u) const override;
+
   // Set or get bias distribution
   void set_bias(std::unique_ptr<PolarAzimuthal> bias)
   {
diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h
index 0d960c33db7..c5823635ebc 100644
--- a/include/openmc/math_functions.h
+++ b/include/openmc/math_functions.h
@@ -201,6 +201,12 @@ std::complex<double> faddeeva(std::complex<double> z);
 //! \return Derivative of Faddeeva function evaluated at z
 std::complex<double> w_derivative(std::complex<double> z, int order);
 
+//! Evaluate relative exponential function
+//!
+//! \param x Real argument
+//! \return (exp(x)-1)/x without loss of precision near 0
+double exprel(double x);
+
 //! Helper function to get index and interpolation function on an incident
 //! energy grid
 //!
diff --git a/include/openmc/particle.h b/include/openmc/particle.h
index 0f37719b94f..6b92b12f190 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -37,6 +37,8 @@ class Particle : public ParticleData {
   //==========================================================================
   // Methods
 
+  double mass() const;
+
   double speed() const;
 
   //! create a secondary particle
@@ -64,10 +66,12 @@ class Particle : public ParticleData {
   //! \param src Source site data
   void from_source(const SourceSite* src);
 
+  void initialize_pseudoparticle(Particle& p, Direction u, double E);
+
   // Coarse-grained particle events
-  void event_calculate_xs();
+  void event_calculate_xs(bool is_pseudo = false);
   void event_advance();
-  void event_cross_surface();
+  void event_cross_surface(bool is_pseudo = false);
   void event_collide();
   void event_revive_from_secondary();
   void event_death();
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/random_lcg.h b/include/openmc/random_lcg.h
index 5aecafed3cb..cc6ccfac47c 100644
--- a/include/openmc/random_lcg.h
+++ b/include/openmc/random_lcg.h
@@ -9,11 +9,12 @@ namespace openmc {
 // Module constants.
 //==============================================================================
 
-constexpr int N_STREAMS {4};
+constexpr int N_STREAMS {5};
 constexpr int STREAM_TRACKING {0};
 constexpr int STREAM_SOURCE {1};
 constexpr int STREAM_URR_PTABLE {2};
 constexpr int STREAM_VOLUME {3};
+constexpr int STREAM_NEXT_EVENT {4};
 constexpr int64_t DEFAULT_SEED {1};
 constexpr uint64_t DEFAULT_STRIDE {152917ULL};
 
diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h
index 4fbbc1b626a..98e5910db9b 100644
--- a/include/openmc/reaction_product.h
+++ b/include/openmc/reaction_product.h
@@ -49,6 +49,11 @@ class ReactionProduct {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
+  AngleEnergy& sample_dist(double E_in, uint64_t* seed) const;
+
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const;
+
   ParticleType particle_;      //!< Particle type
   EmissionMode emission_mode_; //!< Emission mode
   double decay_rate_; //!< Decay rate (for delayed neutron precursors) in [1/s]
diff --git a/include/openmc/secondary_correlated.h b/include/openmc/secondary_correlated.h
index 6905c38e369..e7c933bb614 100644
--- a/include/openmc/secondary_correlated.h
+++ b/include/openmc/secondary_correlated.h
@@ -40,6 +40,9 @@ class CorrelatedAngleEnergy : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  Distribution& sample_dist(double E_in, double& E_out, uint64_t* seed) const;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
   // energy property
   vector<double>& energy() { return energy_; }
diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h
index 83806d35248..0168ee13ced 100644
--- a/include/openmc/secondary_kalbach.h
+++ b/include/openmc/secondary_kalbach.h
@@ -31,6 +31,10 @@ class KalbachMann : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  void sample_params(double E_in, double& E_out, double& km_a, double& km_r,
+    uint64_t* seed) const;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) 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 efb4fd75ba1..dad896cd740 100644
--- a/include/openmc/secondary_nbody.h
+++ b/include/openmc/secondary_nbody.h
@@ -27,6 +27,9 @@ class NBodyPhaseSpace : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy(double E_in, uint64_t* seed) const;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) 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 5b18902afbb..b7de6ab4e61 100644
--- a/include/openmc/secondary_thermal.h
+++ b/include/openmc/secondary_thermal.h
@@ -6,6 +6,7 @@
 
 #include "openmc/angle_energy.h"
 #include "openmc/endf.h"
+#include "openmc/search.h"
 #include "openmc/secondary_correlated.h"
 #include "openmc/vector.h"
 
@@ -32,6 +33,8 @@ class CoherentElasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
@@ -55,6 +58,8 @@ class IncoherentElasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   double debye_waller_;
@@ -80,6 +85,8 @@ class IncoherentElasticAEDiscrete : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const vector<double>& energy_;  //!< Energies at which cosines are tabulated
@@ -106,6 +113,9 @@ class IncoherentInelasticAEDiscrete : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  void sample_params(double E_in, double& E_out, int& j, uint64_t* seed) const;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const vector<double>& energy_; //!< Incident energies
@@ -134,6 +144,10 @@ class IncoherentInelasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  void sample_params(double E_in, double& E_out, double& f, int& l, int& j,
+    uint64_t* seed) const;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   //! Secondary energy/angle distribution
@@ -169,6 +183,9 @@ class MixedElasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  const AngleEnergy& sample_dist(double E_in, uint64_t* seed) const;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   CoherentElasticAE coherent_dist_;         //!< Coherent distribution
@@ -178,6 +195,52 @@ class MixedElasticAE : public AngleEnergy {
   const Function1D& incoherent_xs_;      //!< Polymorphic ref. to incoherent XS
 };
 
+struct DoubleVector {
+  double data;
+  const double& operator[](size_t index) const { return data; }
+};
+
+template<typename T>
+double get_pdf_discrete(const vector<double>& mu, const T& w, double mu_0)
+{
+  // Make sure mu is in range [-1,1]
+  if (std::abs(mu_0) > 1.0)
+    mu_0 = std::copysign(1.0, mu_0);
+  double a0;
+  double a1;
+  double b0;
+  double b1;
+  int32_t ai = -1;
+  int32_t bi = -1;
+  if (mu_0 > mu[0]) {
+    ai = lower_bound_index(mu.begin(), mu.end(), mu_0);
+    a0 = mu[ai];
+    a1 = (ai > 1) ? mu[ai - 1] : -1.0;
+  } else {
+    a0 = -1.0;
+    a1 = -1.0;
+  }
+  if (mu_0 < mu[mu.size() - 1]) {
+    bi = upper_bound_index(mu.begin(), mu.end(), mu_0);
+    b0 = mu[bi];
+    b1 = (bi < mu.size() - 1) ? mu[bi + 1] : 1.0;
+  } else {
+    b0 = 1.0;
+    b1 = 1.0;
+  }
+
+  //  Calculate Delta_a and Delta_b
+  double delta_a = 0.5 * std::min(b0 - a0, a0 - a1);
+  double delta_b = 0.5 * std::min(b1 - b0, b0 - a0);
+
+  if (mu_0 < a0 + delta_a)
+    return w[ai] / (2.0 * delta_a);
+  else if (mu_0 + delta_b < b0)
+    return w[bi] / (2.0 * delta_b);
+  else
+    return 0.0;
+}
+
 } // namespace openmc
 
 #endif // OPENMC_SECONDARY_THERMAL_H
diff --git a/include/openmc/secondary_uncorrelated.h b/include/openmc/secondary_uncorrelated.h
index 3afa3d9ceb7..d971231b277 100644
--- a/include/openmc/secondary_uncorrelated.h
+++ b/include/openmc/secondary_uncorrelated.h
@@ -31,6 +31,8 @@ class UncorrelatedAngleEnergy : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
   // Accessors
   AngleDistribution& angle() { return angle_; }
diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h
index 9a6cf1b2131..42ab43c1fca 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;
@@ -92,6 +94,9 @@ void broadcast_results();
 
 void free_memory_simulation();
 
+//! Simulate a single pseudoparticle history
+double transport_pseudoparticle(Particle& p, double total_distance);
+
 //! Simulate a single particle history (and all generated secondary particles,
 //!  if enabled), from birth to death
 void transport_history_based_single_particle(Particle& p);
diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h
index 65098597a5d..e900d35af9d 100644
--- a/include/openmc/tallies/filter.h
+++ b/include/openmc/tallies/filter.h
@@ -39,6 +39,7 @@ enum class FilterType {
   MUSURFACE,
   PARENT_NUCLIDE,
   PARTICLE,
+  POINT,
   POLAR,
   SPHERICAL_HARMONICS,
   SPATIAL_LEGENDRE,
diff --git a/include/openmc/tallies/filter_point.h b/include/openmc/tallies/filter_point.h
new file mode 100644
index 00000000000..e2f3908d746
--- /dev/null
+++ b/include/openmc/tallies/filter_point.h
@@ -0,0 +1,57 @@
+#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);
+  void reset_indices();
+
+private:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  vector<std::pair<Position, double>> detectors_;
+  vector<vector<int32_t>> indices_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_PARTICLE_H
diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h
index 374daff92a0..dc20505d39d 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"
@@ -12,6 +13,7 @@
 #include "xtensor/xfixed.hpp"
 #include "xtensor/xtensor.hpp"
 
+#include <set>
 #include <string>
 #include <unordered_map>
 
@@ -210,11 +212,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<int> 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 c3ab779e6a1..b3b653383d8 100644
--- a/include/openmc/tallies/tally_scoring.h
+++ b/include/openmc/tallies/tally_scoring.h
@@ -1,12 +1,19 @@
 #ifndef OPENMC_TALLIES_TALLY_SCORING_H
 #define OPENMC_TALLIES_TALLY_SCORING_H
 
+#include "openmc/nuclide.h"
 #include "openmc/particle.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
+#include "openmc/thermal.h"
 
 namespace openmc {
 
+namespace simulation {
+extern thread_local Particle pseudoparticle;
+extern thread_local int i_det;
+} // namespace simulation
+
 //==============================================================================
 //! An iterator over all combinations of a tally's matching filter bins.
 //
@@ -114,6 +121,15 @@ void score_surface_tally(Particle& p, const vector<int>& tallies);
 //! \param tallies A vector of the indices of the tallies to score to
 void score_pulse_height_tally(Particle& p, const vector<int>& tallies);
 
+void score_point_tally(Particle& p, int i_nuclide, const Reaction& rx,
+  int i_product, Direction* v_t);
+
+void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab,
+  const NuclideMicroXS& micro);
+
+void score_point_tally(SourceSite& site, int source_index);
+
+void score_pseudoparticle_tally(Particle& p, double mfp, double pdf);
 } // namespace openmc
 
 #endif // OPENMC_TALLIES_TALLY_SCORING_H
diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h
index de0767d0af0..73c343ad6a7 100644
--- a/include/openmc/thermal.h
+++ b/include/openmc/thermal.h
@@ -51,7 +51,11 @@ class ThermalData {
   //! \param[out] mu Outgoing scattering angle cosine
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(const NuclideMicroXS& micro_xs, double E_in, double* E_out,
-    double* mu, uint64_t* seed);
+    double* mu, uint64_t* seed) const;
+  AngleEnergy& sample_dist(
+    const NuclideMicroXS& micro_xs, double E, uint64_t* seed) const;
+  double sample_energy_and_pdf(const NuclideMicroXS& micro_xs, double E_in,
+    double mu, double& E_out, uint64_t* seed) const;
 
 private:
   struct Reaction {
diff --git a/openmc/filter.py b/openmc/filter.py
index e7c54f2f0ee..572b22270f8 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -22,11 +22,11 @@
 
 _FILTER_TYPES = (
     'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy',
-    'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
-    'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre',
-    'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance',
-    'collision', 'time', 'parentnuclide', 'weight', 'meshborn', 'meshsurface',
-    'meshmaterial',
+    'energyout', 'mu', 'musurface', 'point', 'polar', 'azimuthal', 'distribcell', 
+    'delayedgroup', 'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 
+    'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 
+    'cellinstance', 'collision', 'time', 'parentnuclide', 'weight', 'meshborn', 
+    'meshsurface', 'meshmaterial',
 )
 
 _CURRENT_NAMES = (
@@ -787,6 +787,66 @@ 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
+    
+    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 5cf496107b3..e3cb96d65af 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', 'PolarFilter', 'SphericalHarmonicsFilter',
+    'ParentNuclideFilter', 'ParticleFilter', 'PointFilter', 'PolarFilter', 'SphericalHarmonicsFilter',
     'SpatialLegendreFilter', 'SurfaceFilter', 'TimeFilter', 'UniverseFilter',
     'WeightFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters'
 ]
@@ -603,6 +603,8 @@ def bins(self):
             self._index, particle_i.ctypes.data_as(POINTER(c_int)))
         return [ParticleType(i) for i in particle_i]
 
+class PointFilter(Filter):
+    filter_type = 'point'
 
 class PolarFilter(Filter):
     filter_type = 'polar'
diff --git a/src/chain.cpp b/src/chain.cpp
index a60ef12cd6c..85d9c119791 100644
--- a/src/chain.cpp
+++ b/src/chain.cpp
@@ -74,6 +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
+{
+  E_out = photon_energy_->sample(seed).first;
+  return 0.5;
+}
+
 //==============================================================================
 // Global variables
 //==============================================================================
diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp
index 50f1aca112b..66accac25ae 100644
--- a/src/distribution_angle.cpp
+++ b/src/distribution_angle.cpp
@@ -83,4 +83,19 @@ double AngleDistribution::sample(double E, uint64_t* seed) const
   return mu;
 }
 
+double AngleDistribution::evaluate(double E, double mu) const
+{
+  // Find energy bin and calculate interpolation factor
+  int i;
+  double r;
+  get_energy_index(energy_, E, i, r);
+
+  double pdf = 0.0;
+  if (r > 0.0)
+    pdf += r * distribution_[i + 1]->evaluate(mu);
+  if (r < 1.0)
+    pdf += (1.0 - r) * distribution_[i]->evaluate(mu);
+  return pdf;
+}
+
 } // namespace openmc
diff --git a/src/distribution_multi.cpp b/src/distribution_multi.cpp
index 857e1c30b40..8ceb712ac08 100644
--- a/src/distribution_multi.cpp
+++ b/src/distribution_multi.cpp
@@ -44,6 +44,7 @@ UnitSphereDistribution::UnitSphereDistribution(pugi::xml_node node)
       fatal_error("Angular distribution reference direction must have "
                   "three parameters specified.");
     u_ref_ = Direction(u_ref.data());
+    u_ref_ /= u_ref_.norm();
   }
 }
 
@@ -65,6 +66,7 @@ PolarAzimuthal::PolarAzimuthal(pugi::xml_node node)
       fatal_error("Angular distribution reference v direction must have "
                   "three parameters specified.");
     v_ref_ = Direction(v_ref.data());
+    v_ref_ /= v_ref_.norm();
   }
   w_ref_ = u_ref_.cross(v_ref_);
   if (check_for_node(node, "mu")) {
@@ -116,6 +118,13 @@ std::pair<Direction, double> PolarAzimuthal::sample_impl(
     weight};
 }
 
+double PolarAzimuthal::evaluate(Direction u) const
+{
+  double mu = u.dot(u_ref_);
+  double phi = std::acos(u.dot(v_ref_) / std::sqrt(1 - mu * mu));
+  return mu_->evaluate(mu) * phi_->evaluate(phi);
+}
+
 //==============================================================================
 // Isotropic implementation
 //==============================================================================
@@ -157,6 +166,11 @@ std::pair<Direction, double> Isotropic::sample(uint64_t* seed) const
   }
 }
 
+double Isotropic::evaluate(Direction u) const
+{
+  return 1.0 / (4.0 * PI);
+}
+
 //==============================================================================
 // Monodirectional implementation
 //==============================================================================
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 9473f928b2c..20ca0fb52c5 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -919,6 +919,15 @@ std::complex<double> w_derivative(std::complex<double> z, int order)
   }
 }
 
+double exprel(double x)
+{
+  if (std::abs(x) < 1e-16)
+    return 1.0;
+  else {
+    return std::expm1(x) / x;
+  }
+}
+
 // Helper function to get index and interpolation function on an incident energy
 // grid
 void get_energy_index(
diff --git a/src/particle.cpp b/src/particle.cpp
index b9f9c8f86b6..337914de138 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -43,23 +43,25 @@ namespace openmc {
 // Particle implementation
 //==============================================================================
 
+double Particle::mass() const
+{
+  // Determine mass in eV/c^2
+  switch (this->type()) {
+  case ParticleType::neutron:
+    return MASS_NEUTRON_EV;
+  case ParticleType::photon:
+    return 0.0;
+  case ParticleType::electron:
+  case ParticleType::positron:
+    return MASS_ELECTRON_EV;
+  }
+  UNREACHABLE();
+}
+
 double Particle::speed() const
 {
   if (settings::run_CE) {
-    // Determine mass in eV/c^2
-    double mass;
-    switch (this->type()) {
-    case ParticleType::neutron:
-      mass = MASS_NEUTRON_EV;
-      break;
-    case ParticleType::photon:
-      mass = 0.0;
-      break;
-    case ParticleType::electron:
-    case ParticleType::positron:
-      mass = MASS_ELECTRON_EV;
-      break;
-    }
+    double mass = this->mass();
     // Equivalent to C * sqrt(1-(m/(m+E))^2) without problem at E<<m:
     return C_LIGHT * std::sqrt(this->E() * (this->E() + 2 * mass)) /
            (this->E() + mass);
@@ -157,10 +159,46 @@ void Particle::from_source(const SourceSite* src)
   }
 }
 
-void Particle::event_calculate_xs()
+void Particle::initialize_pseudoparticle(
+  Particle& p, Direction u_new, double E_new)
+{
+  // Reset some attributes
+  clear();
+  surface() = SURFACE_NONE;
+  cell_born() = C_NONE;
+  material() = C_NONE;
+  n_collision() = 0;
+  fission() = false;
+  zero_flux_derivs();
+  lifetime() = 0.0;
+
+  // Copy attributes from particle
+  type() = p.type();
+  wgt() = p.wgt();
+  wgt_last() = p.wgt();
+  r() = p.r();
+  u() = u_new;
+  r_born() = p.r();
+  r_last_current() = p.r();
+  r_last() = p.r();
+  u_last() = p.u();
+  if (settings::run_CE) {
+    E() = E_new;
+    g() = 0;
+  } else {
+    g() = static_cast<int>(E_new);
+    g_last() = static_cast<int>(E_new);
+    E() = data::mg.energy_bin_avg_[g()];
+  }
+  E_last() = E();
+  time() = p.time();
+  time_last() = p.time();
+}
+
+void Particle::event_calculate_xs(bool is_pseudo)
 {
   // Set the random number stream
-  stream() = STREAM_TRACKING;
+  stream() = is_pseudo ? STREAM_NEXT_EVENT : STREAM_TRACKING;
 
   // Store pre-collision particle properties
   wgt_last() = wgt();
@@ -196,7 +234,7 @@ void Particle::event_calculate_xs()
   }
 
   // Write particle track.
-  if (write_track())
+  if (write_track() && !is_pseudo)
     write_particle_track(*this);
 
   if (settings::check_overlaps)
@@ -285,7 +323,7 @@ void Particle::event_advance()
   }
 }
 
-void Particle::event_cross_surface()
+void Particle::event_cross_surface(bool is_pseudo)
 {
   // Saving previous cell data
   for (int j = 0; j < n_coord(); ++j) {
@@ -309,21 +347,21 @@ void Particle::event_cross_surface()
     // Particle crosses surface
     const auto& surf {model::surfaces[surface_index()].get()};
     // If BC, add particle to surface source before crossing surface
-    if (surf->surf_source_ && surf->bc_) {
+    if (surf->surf_source_ && surf->bc_ && !is_pseudo) {
       add_surf_source_to_bank(*this, *surf);
     }
     this->cross_surface(*surf);
     // If no BC, add particle to surface source after crossing surface
-    if (surf->surf_source_ && !surf->bc_) {
+    if (surf->surf_source_ && !surf->bc_ && !is_pseudo) {
       add_surf_source_to_bank(*this, *surf);
     }
-    if (settings::weight_window_checkpoint_surface) {
+    if (settings::weight_window_checkpoint_surface && !is_pseudo) {
       apply_weight_windows(*this);
     }
     event() = TallyEvent::SURFACE;
   }
   // Score cell to cell partial currents
-  if (!model::active_surface_tallies.empty()) {
+  if (!model::active_surface_tallies.empty() && !is_pseudo) {
     score_surface_tally(*this, model::active_surface_tallies);
   }
 }
diff --git a/src/physics.cpp b/src/physics.cpp
index 41509af97be..e62029e6f17 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -25,6 +25,7 @@
 #include "openmc/simulation.h"
 #include "openmc/string_utils.h"
 #include "openmc/tallies/tally.h"
+#include "openmc/tallies/tally_scoring.h"
 #include "openmc/thermal.h"
 #include "openmc/weight_windows.h"
 
@@ -732,7 +733,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_;
   }
 
@@ -778,6 +779,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(p, i_nuclide, rx, 0, &v_t);
+  }
+
   // Sample scattering angle, checking if angle distribution is present (assume
   // isotropic otherwise)
   double mu_cm;
@@ -825,8 +830,12 @@ 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];
+  sab.sample(micro, p.E(), &E_out, &p.mu(), p.current_seed());
+
+  if (!model::active_point_tallies.empty()) {
+    score_point_tally(p, i_nuclide, sab, micro);
+  }
 
   // Set energy to outgoing, change direction of particle
   p.E() = E_out;
@@ -1091,6 +1100,10 @@ void sample_fission_neutron(
     site->delayed_group = 0;
   }
 
+  if (!model::active_point_tallies.empty()) {
+    score_point_tally(p, i_nuclide, rx, site->delayed_group, nullptr);
+  }
+
   // sample from prompt neutron energy distribution
   int n_sample = 0;
   double mu;
@@ -1116,8 +1129,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();
 
@@ -1126,13 +1142,17 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
   double mu;
   rx.products_[0].sample(E_in, E, mu, p.current_seed());
 
+  if (!model::active_point_tallies.empty()) {
+    score_point_tally(p, i_nuclide, rx, 0, nullptr);
+  }
+
   // 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));
 
diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp
index 3ba2c0cfb05..6cfd164a01d 100644
--- a/src/reaction_product.cpp
+++ b/src/reaction_product.cpp
@@ -106,8 +106,7 @@ ReactionProduct::ReactionProduct(const ChainNuclide::Product& product)
     make_unique<DecayPhotonAngleEnergy>(chain_nuc->photon_energy()));
 }
 
-void ReactionProduct::sample(
-  double E_in, double& E_out, double& mu, uint64_t* seed) const
+AngleEnergy& ReactionProduct::sample_dist(double E_in, uint64_t* seed) const
 {
   auto n = applicability_.size();
   if (n > 1) {
@@ -118,15 +117,25 @@ void ReactionProduct::sample(
       prob += applicability_[i](E_in);
 
       // If i-th distribution is sampled, sample energy from the distribution
-      if (c <= prob) {
-        distribution_[i]->sample(E_in, E_out, mu, seed);
-        break;
-      }
+      if (c <= prob)
+        return *distribution_[i];
     }
   } else {
     // If only one distribution is present, go ahead and sample it
-    distribution_[0]->sample(E_in, E_out, mu, seed);
+    return *distribution_[0];
   }
 }
 
+void ReactionProduct::sample(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const
+{
+  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
+{
+  return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed);
+}
+
 } // namespace openmc
diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp
index e820419b484..dfde3f1e03b 100644
--- a/src/secondary_correlated.cpp
+++ b/src/secondary_correlated.cpp
@@ -153,9 +153,8 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group)
     distribution_.push_back(std::move(d));
   } // incoming energies
 }
-
-void CorrelatedAngleEnergy::sample(
-  double E_in, double& E_out, double& mu, uint64_t* seed) const
+Distribution& CorrelatedAngleEnergy::sample_dist(
+  double E_in, double& E_out, uint64_t* seed) const
 {
   // Find energy bin and calculate interpolation factor
   int i;
@@ -247,10 +246,22 @@ void CorrelatedAngleEnergy::sample(
   // Find correlated angular distribution for closest outgoing energy bin
   if (r1 - c_k < c_k1 - r1 ||
       distribution_[l].interpolation == Interpolation::histogram) {
-    mu = distribution_[l].angle[k]->sample(seed).first;
+    return *distribution_[l].angle[k];
   } else {
-    mu = distribution_[l].angle[k + 1]->sample(seed).first;
+    return *distribution_[l].angle[k + 1];
   }
 }
 
+void CorrelatedAngleEnergy::sample(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const
+{
+  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
+{
+  return sample_dist(E_in, E_out, seed).evaluate(mu);
+}
+
 } // namespace openmc
diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp
index 6ac91e665b9..43a40c0cd4e 100644
--- a/src/secondary_kalbach.cpp
+++ b/src/secondary_kalbach.cpp
@@ -115,8 +115,8 @@ KalbachMann::KalbachMann(hid_t group)
   } // incoming energies
 }
 
-void KalbachMann::sample(
-  double E_in, double& E_out, double& mu, uint64_t* seed) const
+void KalbachMann::sample_params(
+  double E_in, double& E_out, double& km_a, double& km_r, uint64_t* seed) const
 {
   // Find energy bin and calculate interpolation factor
   int i;
@@ -171,7 +171,6 @@ void KalbachMann::sample(
 
   double E_l_k = distribution_[l].e_out[k];
   double p_l_k = distribution_[l].p[k];
-  double km_r, km_a;
   if (distribution_[l].interpolation == Interpolation::histogram) {
     // Histogram interpolation
     if (p_l_k > 0.0 && k >= n_discrete) {
@@ -217,6 +216,13 @@ void KalbachMann::sample(
       E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
     }
   }
+}
+
+void KalbachMann::sample(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const
+{
+  double km_r, km_a;
+  sample_params(E_in, E_out, km_a, km_r, seed);
 
   // Sampled correlated angle from Kalbach-Mann parameters
   if (prn(seed) > km_r) {
@@ -227,5 +233,15 @@ 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 km_r, km_a;
+  sample_params(E_in, E_out, km_a, km_r, seed);
+
+  // https://docs.openmc.org/en/latest/methods/neutron_physics.html#equation-KM-pdf-angle
+  return km_a / (2 * std::sinh(km_a)) *
+         (std::cosh(km_a * mu) + km_r * std::sinh(km_a * mu));
+}
 
 } // namespace openmc
diff --git a/src/secondary_nbody.cpp b/src/secondary_nbody.cpp
index da0bb81c471..72f0b0b92d0 100644
--- a/src/secondary_nbody.cpp
+++ b/src/secondary_nbody.cpp
@@ -22,13 +22,8 @@ NBodyPhaseSpace::NBodyPhaseSpace(hid_t group)
   read_attribute(group, "q_value", Q_);
 }
 
-void NBodyPhaseSpace::sample(
-  double E_in, double& E_out, double& mu, uint64_t* seed) const
+double NBodyPhaseSpace::sample_energy(double E_in, uint64_t* seed) const
 {
-  // By definition, the distribution of the angle is isotropic for an N-body
-  // phase space distribution
-  mu = uniform_distribution(-1., 1., seed);
-
   // Determine E_max parameter
   double Ap = mass_ratio_;
   double E_max = (Ap - 1.0) / Ap * (A_ / (A_ + 1.0) * E_in + Q_);
@@ -59,12 +54,29 @@ void NBodyPhaseSpace::sample(
         std::log(r5) * std::pow(std::cos(PI / 2.0 * r6), 2);
     break;
   default:
-    throw std::runtime_error {"N-body phase space with >5 bodies."};
+    fatal_error("N-body phase space with >5 bodies.");
   }
 
   // Now determine v and E_out
   double v = x / (x + y);
-  E_out = E_max * v;
+  return E_max * v;
+}
+
+void NBodyPhaseSpace::sample(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const
+{
+  // By definition, the distribution of the angle is isotropic for an N-body
+  // phase space distribution
+  mu = uniform_distribution(-1., 1., seed);
+
+  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
+{
+  E_out = sample_energy(E_in, seed);
+  return 0.5;
 }
 
 } // namespace openmc
diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp
index 030d398aabe..240bbf8cf36 100644
--- a/src/secondary_thermal.cpp
+++ b/src/secondary_thermal.cpp
@@ -4,6 +4,7 @@
 #include "openmc/math_functions.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
+#include "openmc/vector.h"
 
 #include "xtensor/xview.hpp"
 
@@ -12,6 +13,12 @@
 
 namespace openmc {
 
+double get_pdf_discrete(const vector<double>& mu, double mu_0)
+{
+  DoubleVector w {1.0 / mu.size()};
+  return get_pdf_discrete(mu, w, mu_0);
+}
+
 //==============================================================================
 // CoherentElasticAE implementation
 //==============================================================================
@@ -42,6 +49,39 @@ 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
+{
+  // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
+
+  double pdf;
+  E_out = E_in;
+  const auto& energies {xs_.bragg_edges()};
+  const auto& factors = xs_.factors();
+
+  if (E_in < energies.front() || E_in > energies.back()) {
+    return 0;
+  }
+
+  const int n = upper_bound_index(energies.begin(), energies.end(), E_in);
+
+  vector<double> mu_vector;
+  mu_vector.reserve(n);
+
+  std::transform(energies.rbegin() + n - 1, energies.rend(),
+    std::back_inserter(mu_vector),
+    [E_in](double Ei) { return 1 - 2 * Ei / E_in; });
+
+  vector<double> weights;
+  weights.reserve(n);
+
+  weights.emplace_back(factors[0] / factors[n]);
+  for (int i = 1; i <= n; ++i) {
+    weights.emplace_back((factors[i] - factors[i - 1]) / factors[n]);
+  }
+  return get_pdf_discrete(mu_vector, weights, mu);
+}
+
 //==============================================================================
 // IncoherentElasticAE implementation
 //==============================================================================
@@ -54,12 +94,21 @@ IncoherentElasticAE::IncoherentElasticAE(hid_t group)
 void IncoherentElasticAE::sample(
   double E_in, double& E_out, double& mu, uint64_t* seed) const
 {
+  E_out = E_in;
+
   // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
   double c = 2 * E_in * debye_waller_;
   mu = std::log(1.0 + prn(seed) * (std::exp(2.0 * c) - 1)) / c - 1.0;
-
-  // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4)
+}
+double IncoherentElasticAE::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
   E_out = E_in;
+
+  // 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));
 }
 
 //==============================================================================
@@ -116,6 +165,28 @@ 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
+{
+  // Get index and interpolation factor for elastic grid
+  int i;
+  double f;
+  get_energy_index(energy_, E_in, i, f);
+  // Energy doesn't change in elastic scattering
+  E_out = E_in;
+  int n_mu = mu_out_.shape()[1];
+
+  std::vector<double> mu_vector;
+  mu_vector.reserve(n_mu);
+
+  for (int k = 0; k < n_mu; ++k) {
+    double mu_k = mu_out_(i, k) + f * (mu_out_(i + 1, k) - mu_out_(i, k));
+    mu_vector.push_back(mu_k);
+  }
+
+  return get_pdf_discrete(mu_vector, mu);
+}
+
 //==============================================================================
 // IncoherentInelasticAEDiscrete implementation
 //==============================================================================
@@ -129,8 +200,8 @@ IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete(
   read_dataset(group, "skewed", skewed_);
 }
 
-void IncoherentInelasticAEDiscrete::sample(
-  double E_in, double& E_out, double& mu, uint64_t* seed) const
+void IncoherentInelasticAEDiscrete::sample_params(
+  double E_in, double& E_out, int& j, uint64_t* seed) const
 {
   // Get index and interpolation factor for inelastic grid
   int i;
@@ -144,7 +215,6 @@ void IncoherentInelasticAEDiscrete::sample(
   // for the second and second to last bins, relative to a normal bin
   // probability of 1). Otherwise, each bin is equally probable.
 
-  int j;
   int n = energy_out_.shape()[1];
   if (!skewed_) {
     // All bins equally likely
@@ -176,6 +246,18 @@ void IncoherentInelasticAEDiscrete::sample(
 
   // Outgoing energy
   E_out = (1 - f) * E_ij + f * E_i1j;
+}
+
+void IncoherentInelasticAEDiscrete::sample(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const
+{
+  // Get index and interpolation factor for inelastic grid
+  int i;
+  double f;
+  get_energy_index(energy_, E_in, i, f);
+
+  int j;
+  sample_params(E_in, E_out, j, seed);
 
   // Sample outgoing cosine bin
   int m = mu_out_.shape()[2];
@@ -189,6 +271,30 @@ 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
+{
+  // Get index and interpolation factor for inelastic grid
+  int i;
+  double f;
+  get_energy_index(energy_, E_in, i, f);
+  int j;
+  sample_params(E_in, E_out, j, seed);
+
+  int m = mu_out_.shape()[2];
+  std::vector<double> mu_vector;
+  mu_vector.reserve(m);
+
+  for (int k = 0; k < m; ++k) {
+    double mu_ijk = mu_out_(i, j, k);
+    double mu_i1jk = mu_out_(i + 1, j, k);
+    double mu_k = (1 - f) * mu_ijk + f * mu_i1jk;
+    mu_vector.push_back(mu_k);
+  }
+
+  return get_pdf_discrete(mu_vector, mu);
+}
+
 //==============================================================================
 // IncoherentInelasticAE implementation
 //==============================================================================
@@ -231,24 +337,23 @@ IncoherentInelasticAE::IncoherentInelasticAE(hid_t group)
   }
 }
 
-void IncoherentInelasticAE::sample(
-  double E_in, double& E_out, double& mu, uint64_t* seed) const
+void IncoherentInelasticAE::sample_params(
+  double E_in, double& E_out, double& f, int& l, int& j, uint64_t* seed) const
 {
   // Get index and interpolation factor for inelastic grid
   int i;
-  double f;
-  get_energy_index(energy_, E_in, i, f);
+  double f0;
+  get_energy_index(energy_, E_in, i, f0);
 
   // Pick closer energy based on interpolation factor
-  int l = f > 0.5 ? i + 1 : i;
+  l = f0 > 0.5 ? i + 1 : i;
 
   // Determine outgoing energy bin
   // (First reset n_energy_out to the right value)
-  auto n = distribution_[l].n_e_out;
+  int n = distribution_[l].n_e_out;
   double r1 = prn(seed);
   double c_j = distribution_[l].e_out_cdf[0];
   double c_j1;
-  std::size_t j;
   for (j = 0; j < n - 1; ++j) {
     c_j1 = distribution_[l].e_out_cdf[j + 1];
     if (r1 < c_j1)
@@ -286,6 +391,15 @@ void IncoherentInelasticAE::sample(
     E_out += E_in - E_l;
   }
 
+  f = (r1 - c_j) / (c_j1 - c_j);
+}
+void IncoherentInelasticAE::sample(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const
+{
+  double f;
+  int l, j;
+  sample_params(E_in, E_out, f, l, j, seed);
+
   // Sample outgoing cosine bin
   int n_mu = distribution_[l].mu.shape()[1];
   std::size_t k = prn(seed) * n_mu;
@@ -294,7 +408,6 @@ void IncoherentInelasticAE::sample(
   // a bin of width 0.5*min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
   // discrete mu value itself.
   const auto& mu_l = distribution_[l].mu;
-  f = (r1 - c_j) / (c_j1 - c_j);
 
   // Interpolate kth mu value between distributions at energies j and j+1
   mu = mu_l(j, k) + f * (mu_l(j + 1, k) - mu_l(j, k));
@@ -318,6 +431,25 @@ 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 f;
+  int l, j;
+  sample_params(E_in, E_out, f, l, j, seed);
+
+  int n_mu = distribution_[l].mu.shape()[1];
+  const auto& mu_l = distribution_[l].mu;
+  std::vector<double> mu_vector;
+
+  for (int k = 0; k < n_mu; ++k) {
+    double mu_k = mu_l(j, k) + f * (mu_l(j + 1, k) - mu_l(j, k));
+    mu_vector.push_back(mu_k);
+  }
+
+  return get_pdf_discrete(mu_vector, mu);
+}
+
 //==============================================================================
 // MixedElasticAE implementation
 //==============================================================================
@@ -340,18 +472,30 @@ MixedElasticAE::MixedElasticAE(
   close_group(incoherent_group);
 }
 
-void MixedElasticAE::sample(
-  double E_in, double& E_out, double& mu, uint64_t* seed) const
+const AngleEnergy& MixedElasticAE::sample_dist(
+  double E_in, uint64_t* seed) const
 {
   // Evaluate coherent and incoherent elastic cross sections
   double xs_coh = coherent_xs_(E_in);
   double xs_incoh = incoherent_xs_(E_in);
 
   if (prn(seed) * (xs_coh + xs_incoh) < xs_coh) {
-    coherent_dist_.sample(E_in, E_out, mu, seed);
+    return coherent_dist_;
   } else {
-    incoherent_dist_->sample(E_in, E_out, mu, seed);
+    return *incoherent_dist_;
   }
 }
 
+void MixedElasticAE::sample(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const
+{
+  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
+{
+  return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed);
+}
+
 } // namespace openmc
diff --git a/src/secondary_uncorrelated.cpp b/src/secondary_uncorrelated.cpp
index 5cbb76fb9b9..ec2af710258 100644
--- a/src/secondary_uncorrelated.cpp
+++ b/src/secondary_uncorrelated.cpp
@@ -65,4 +65,22 @@ 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
+{
+  // Sample outgoing energy
+  if (energy_ != nullptr) {
+    E_out = energy_->sample(E_in, seed);
+  } else {
+    E_out = E_in;
+  }
+
+  if (!angle_.empty()) {
+    return angle_.evaluate(E_in, mu);
+  } else {
+    // no angle distribution given => assume isotropic for all energies
+    return 0.5;
+  }
+}
+
 } // namespace openmc
diff --git a/src/simulation.cpp b/src/simulation.cpp
index b536ae5881f..453bb6f20dc 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -7,10 +7,12 @@
 #include "openmc/eigenvalue.h"
 #include "openmc/error.h"
 #include "openmc/event.h"
+#include "openmc/geometry.h"
 #include "openmc/geometry_aux.h"
 #include "openmc/ifp.h"
 #include "openmc/material.h"
 #include "openmc/message_passing.h"
+#include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/output.h"
 #include "openmc/particle.h"
@@ -22,6 +24,7 @@
 #include "openmc/tallies/derivative.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
+#include "openmc/tallies/tally_scoring.h"
 #include "openmc/tallies/trigger.h"
 #include "openmc/timer.h"
 #include "openmc/track_output.h"
@@ -310,6 +313,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;
@@ -820,6 +824,43 @@ void transport_history_based_single_particle(Particle& p)
   p.event_death();
 }
 
+double transport_pseudoparticle(Particle& p, double total_distance)
+{
+  double time_cutoff = settings::time_cutoff[static_cast<int>(p.type())];
+  double speed = p.speed();
+  p.event_calculate_xs(true);
+
+  double mfp = 0.0;
+  while (total_distance > 0.0) {
+    p.boundary() = distance_to_boundary(p);
+
+    double distance_cutoff =
+      (time_cutoff < INFTY) ? (time_cutoff - p.time()) * speed : INFTY;
+
+    double distance =
+      std::min({p.boundary().distance(), distance_cutoff, total_distance});
+    if (distance == distance_cutoff) {
+      p.wgt() = 0.0;
+      return mfp;
+    }
+
+    // Advance particle in space and time
+    p.move_distance(distance);
+    p.time() += distance / speed;
+    p.lifetime() += distance / speed;
+    total_distance -= distance;
+    mfp += distance * p.macro_xs().total;
+
+    if (distance == p.boundary().distance()) {
+      p.event_cross_surface(true);
+      p.event_calculate_xs(true);
+      if (!settings::run_CE)
+        speed = p.speed();
+    }
+  }
+  return mfp;
+}
+
 void transport_history_based()
 {
 #pragma omp parallel for schedule(runtime)
diff --git a/src/source.cpp b/src/source.cpp
index b30b964d39c..a1067a2bc3f 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -33,6 +33,7 @@
 #include "openmc/simulation.h"
 #include "openmc/state_point.h"
 #include "openmc/string_utils.h"
+#include "openmc/tallies/tally_scoring.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -668,6 +669,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(site, i);
+  }
+
   return site;
 }
 
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 79817981db0..29737b7b23a 100644
--- a/src/tallies/filter.cpp
+++ b/src/tallies/filter.cpp
@@ -31,6 +31,7 @@
 #include "openmc/tallies/filter_musurface.h"
 #include "openmc/tallies/filter_parent_nuclide.h"
 #include "openmc/tallies/filter_particle.h"
+#include "openmc/tallies/filter_point.h"
 #include "openmc/tallies/filter_polar.h"
 #include "openmc/tallies/filter_sph_harm.h"
 #include "openmc/tallies/filter_sptl_legendre.h"
@@ -146,6 +147,8 @@ Filter* Filter::create(const std::string& type, int32_t id)
     return Filter::create<ParentNuclideFilter>(id);
   } else if (type == "particle") {
     return Filter::create<ParticleFilter>(id);
+  } else if (type == "point") {
+    return Filter::create<PointFilter>(id);
   } else if (type == "polar") {
     return Filter::create<PolarFilter>(id);
   } else if (type == "surface") {
diff --git a/src/tallies/filter_point.cpp b/src/tallies/filter_point.cpp
new file mode 100644
index 00000000000..0a9587f6650
--- /dev/null
+++ b/src/tallies/filter_point.cpp
@@ -0,0 +1,92 @@
+#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::reset_indices()
+{
+  indices_.clear();
+  for (auto det : model::active_point_detectors) {
+    vector<int32_t> temp;
+    for (auto i = 0; i < detectors_.size(); i++) {
+      auto [d, r] = detectors_[i];
+      if (d == det)
+        temp.push_back(i);
+    }
+    indices_.push_back(temp);
+  }
+}
+
+void PointFilter::get_all_bins(
+  const Particle& p, TallyEstimator estimator, FilterMatch& match) const
+{
+  double attenuation = p.wgt() / p.wgt_last();
+  for (auto i : indices_[simulation::i_det]) {
+    auto [pos, r] = detectors_[i];
+    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);
+    }
+  }
+}
+
+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/tally.cpp b/src/tallies/tally.cpp
index 6eef1da9cfd..c2f81dd469b 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -11,6 +11,7 @@
 #include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/particle.h"
+#include "openmc/position.h"
 #include "openmc/reaction.h"
 #include "openmc/reaction_product.h"
 #include "openmc/settings.h"
@@ -30,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"
@@ -61,11 +63,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<int> pulse_height_cells;
 vector<double> time_grid;
+std::set<Position> active_point_detectors;
 } // namespace model
 
 namespace simulation {
@@ -535,6 +539,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 surface_present = false;
   bool meshsurface_present = false;
   bool non_cell_energy_present = false;
@@ -551,6 +556,10 @@ void Tally::set_scores(const vector<std::string>& scores)
       cellfrom_present = true;
     } else if (filt->type() == FilterType::CELL) {
       cell_present = true;
+    } else if (filt->type() == FilterType::POINT) {
+      point_present = true;
+      type_ = TallyType::NEXT_EVENT;
+      estimator_ = TallyEstimator::POINT;
     } else if (filt->type() == FilterType::SURFACE) {
       surface_present = true;
     } else if (filt->type() == FilterType::MESH_SURFACE) {
@@ -558,6 +567,19 @@ void Tally::set_scores(const vector<std::string>& scores)
     }
   }
 
+  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
@@ -601,6 +623,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)
@@ -609,6 +634,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 (surface_present || cell_present || cellfrom_present) {
         if (meshsurface_present)
@@ -624,11 +651,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");
@@ -655,6 +686,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;
     }
@@ -1136,6 +1169,14 @@ void add_to_time_grid(vector<double> grid)
   model::time_grid.swap(merged);
 }
 
+//! Add new points to global point_detectors
+//
+//! \param detector Position of new point detector to add
+void add_point_detector(Position& detector)
+{
+  model::active_point_detectors.insert(detector);
+}
+
 void setup_active_tallies()
 {
   model::active_tallies.clear();
@@ -1143,10 +1184,14 @@ void setup_active_tallies()
   model::active_tracklength_tallies.clear();
   model::active_timed_tracklength_tallies.clear();
   model::active_collision_tallies.clear();
+  model::active_point_tallies.clear();
   model::active_meshsurf_tallies.clear();
   model::active_surface_tallies.clear();
   model::active_pulse_height_tallies.clear();
   model::time_grid.clear();
+  model::active_point_detectors.clear();
+
+  std::set<int32_t> particle_filter_ids;
 
   for (auto i = 0; i < model::tallies.size(); ++i) {
     const auto& tally {*model::tallies[i]};
@@ -1187,9 +1232,25 @@ void setup_active_tallies()
       case TallyType::PULSE_HEIGHT:
         model::active_pulse_height_tallies.push_back(i);
         break;
+
+      case TallyType::NEXT_EVENT:
+        switch (tally.estimator_) {
+        case TallyEstimator::POINT:
+          model::active_point_tallies.push_back(i);
+          auto pf = tally.get_filter<PointFilter>();
+          for (auto [det, r] : pf->detectors()) {
+            add_point_detector(det);
+          }
+          particle_filter_ids.insert(model::tally_map[pf->id()]);
+          break;
+        }
       }
     }
   }
+  for (auto id : particle_filter_ids) {
+    auto pf = dynamic_cast<PointFilter*>(model::tally_filters.at(id).get());
+    pf->reset_indices();
+  }
 }
 
 void free_memory_tally()
@@ -1207,10 +1268,12 @@ void free_memory_tally()
   model::active_tracklength_tallies.clear();
   model::active_timed_tracklength_tallies.clear();
   model::active_collision_tallies.clear();
+  model::active_point_tallies.clear();
   model::active_meshsurf_tallies.clear();
   model::active_surface_tallies.clear();
   model::active_pulse_height_tallies.clear();
   model::time_grid.clear();
+  model::active_point_detectors.clear();
 
   model::tally_map.clear();
 }
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index 67e851644a9..4db3bec70f8 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -13,6 +13,7 @@
 #include "openmc/search.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
+#include "openmc/source.h"
 #include "openmc/string_utils.h"
 #include "openmc/tallies/derivative.h"
 #include "openmc/tallies/filter.h"
@@ -24,6 +25,11 @@
 
 namespace openmc {
 
+namespace simulation {
+thread_local Particle pseudoparticle;
+thread_local int i_det;
+} // namespace simulation
+
 //==============================================================================
 // FilterBinIter implementation
 //==============================================================================
@@ -2738,4 +2744,137 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
     p.E_last() = orig_E_last;
   }
 }
+
+void score_pseudoparticle_tally(Particle& p, double mfp, double pdf)
+{
+  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);
+}
+
+void score_point_tally(
+  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_;
+  double E_in = p.E();
+  double vel = std::sqrt(E_in);
+  Direction v_n = vel * p.u();
+  Direction v_cm = v_n / (awr + 1.0);
+  if (v_t != nullptr)
+    v_cm += awr / (awr + 1.0) * (*v_t);
+  Direction u_cm = (v_n - v_cm);
+  u_cm /= u_cm.norm();
+  double mfp = 0;
+
+  auto old_stream = p.stream();
+  p.stream() = STREAM_NEXT_EVENT;
+
+  simulation::i_det = -1;
+  for (auto& det : model::active_point_detectors) {
+    ++simulation::i_det;
+
+    auto u = (det - p.r());
+    double distance = u.norm();
+    u /= distance;
+
+    double mu;
+    if (rx.scatter_in_cm_) {
+      mu = u.dot(u_cm);
+    } else {
+      mu = u.dot(p.u());
+    }
+    double E_out;
+    double pdf = rx.products_[i_product].sample_energy_and_pdf(
+      p.E(), mu, E_out, p.current_seed());
+    if (rx.scatter_in_cm_) {
+      double E_cm = E_out;
+      Direction v_out = std::sqrt(E_cm) * u_cm;
+      v_out += v_cm;
+      E_out = std::pow(v_out.norm(), 2);
+      pdf *= std::sqrt(E_out / E_cm) /
+             (1 - mu / (awr + 1) * std::sqrt(E_in / E_out));
+    }
+    simulation::pseudoparticle.initialize_pseudoparticle(p, u, E_out);
+    double mfp = transport_pseudoparticle(simulation::pseudoparticle, distance);
+    if (!simulation::pseudoparticle.alive())
+      continue;
+    score_pseudoparticle_tally(simulation::pseudoparticle, mfp, pdf);
+  }
+  p.stream() = old_stream;
+}
+
+void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab,
+  const NuclideMicroXS& micro)
+{
+  const auto& nuc {data::nuclides[i_nuclide]};
+  double awr = nuc->awr_;
+  double E_in = p.E();
+  double vel = std::sqrt(E_in);
+  Direction v_n = vel * p.u();
+  Direction v_cm = v_n / (awr + 1.0);
+  Direction u_cm = (v_n - v_cm);
+  u_cm /= u_cm.norm();
+  double mfp = 0;
+
+  auto old_stream = p.stream();
+  p.stream() = STREAM_NEXT_EVENT;
+
+  simulation::i_det = -1;
+  for (auto& det : model::active_point_detectors) {
+    ++simulation::i_det;
+
+    auto u = (det - p.r());
+    double distance = u.norm();
+    u /= distance;
+
+    double mu = u.dot(p.u());
+    double E_out;
+    double pdf =
+      sab.sample_energy_and_pdf(micro, p.E(), mu, E_out, p.current_seed());
+    simulation::pseudoparticle.initialize_pseudoparticle(p, u, E_out);
+    double mfp = transport_pseudoparticle(simulation::pseudoparticle, distance);
+    if (!simulation::pseudoparticle.alive())
+      continue;
+    score_pseudoparticle_tally(simulation::pseudoparticle, mfp, pdf);
+  }
+  p.stream() = old_stream;
+}
+
+void score_point_tally(SourceSite& site, int source_index)
+{
+  double E_out = site.E;
+  Position r = site.r;
+  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 angle = src->angle();
+
+  double mfp = 0;
+
+  simulation::i_det = -1;
+  for (auto& det : model::active_point_detectors) {
+    ++simulation::i_det;
+
+    auto u = (det - r);
+    double distance = u.norm();
+    u /= distance;
+
+    double pdf = angle->evaluate(u);
+    simulation::pseudoparticle.from_source(&site);
+    simulation::pseudoparticle.u() = u;
+    double mfp = transport_pseudoparticle(simulation::pseudoparticle, distance);
+    if (!simulation::pseudoparticle.alive())
+      continue;
+    score_pseudoparticle_tally(simulation::pseudoparticle, mfp, pdf);
+  }
+}
+
 } // namespace openmc
diff --git a/src/thermal.cpp b/src/thermal.cpp
index cbe0983ed65..1a06d1d19a3 100644
--- a/src/thermal.cpp
+++ b/src/thermal.cpp
@@ -296,16 +296,21 @@ void ThermalData::calculate_xs(
   *inelastic = (*inelastic_.xs)(E);
 }
 
-void ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
-  double* E_out, double* mu, uint64_t* seed)
+AngleEnergy& ThermalData::sample_dist(
+  const NuclideMicroXS& micro_xs, double E, uint64_t* seed) const
 {
   // Determine whether inelastic or elastic scattering will occur
   if (prn(seed) < micro_xs.thermal_elastic / micro_xs.thermal) {
-    elastic_.distribution->sample(E, *E_out, *mu, seed);
+    return *elastic_.distribution;
   } else {
-    inelastic_.distribution->sample(E, *E_out, *mu, seed);
+    return *inelastic_.distribution;
   }
+}
 
+void ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
+  double* E_out, double* mu, uint64_t* seed) const
+{
+  sample_dist(micro_xs, E, seed).sample(E, *E_out, *mu, seed);
   // Because of floating-point roundoff, it may be possible for mu to be
   // outside of the range [-1,1). In these cases, we just set mu to exactly
   // -1 or 1
@@ -313,6 +318,13 @@ void ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
     *mu = std::copysign(1.0, *mu);
 }
 
+double ThermalData::sample_energy_and_pdf(const NuclideMicroXS& micro_xs,
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  return sample_dist(micro_xs, E_in, seed)
+    .sample_energy_and_pdf(E_in, mu, E_out, seed);
+}
+
 void free_memory_thermal()
 {
   data::thermal_scatt.clear();