diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90
index 29fd2bf1e0..7c037fcebe 100644
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@@ -51,6 +51,11 @@ module tally
       type(Particle), intent(in) :: p
     end subroutine
 
+    subroutine score_analog_tally_mg(p) bind(C)
+      import Particle
+      type(Particle), intent(in) :: p
+    end subroutine
+
     subroutine score_tracklength_tally(p, distance) bind(C)
       import Particle, C_DOUBLE
       type(Particle) :: p
@@ -2077,141 +2082,6 @@ contains
     end associate
   end subroutine score_all_nuclides
 
-!===============================================================================
-! SCORE_ANALOG_TALLY keeps track of how many events occur in a specified cell,
-! energy range, etc. Note that since these are "analog" tallies, they are only
-! triggered at every collision, not every event
-!===============================================================================
-
-  subroutine score_analog_tally_mg(p)
-
-    type(Particle), intent(in) :: p
-
-    integer :: i, j
-    integer :: i_tally
-    integer :: i_filt
-    integer :: i_bin
-    integer :: k                    ! loop index for nuclide bins
-    integer :: filter_index         ! single index for single bin
-    integer :: i_nuclide            ! index in nuclides array
-    real(8) :: filter_weight        ! combined weight of all filters
-    real(8) :: atom_density
-    logical :: finished             ! found all valid bin combinations
-    type(Material),    pointer :: mat
-
-    ! A loop over all tallies is necessary because we need to simultaneously
-    ! determine different filter bins for the same tally in order to score to it
-
-    TALLY_LOOP: do i = 1, active_analog_tallies_size()
-      ! Get index of tally and pointer to tally
-      i_tally = active_analog_tallies_data(i)
-      associate (t => tallies(i_tally) % obj)
-
-      ! Find all valid bins in each filter if they have not already been found
-      ! for a previous tally.
-      do j = 1, t % n_filters()
-        i_filt = t % filter(j)
-        if (.not. filter_matches(i_filt) % bins_present) then
-          call filter_matches(i_filt) % bins_clear()
-          call filter_matches(i_filt) % weights_clear()
-          call filters(i_filt) % obj % get_all_bins(p, t % estimator(), &
-               filter_matches(i_filt))
-          filter_matches(i_filt) % bins_present = .true.
-        end if
-        ! If there are no valid bins for this filter, then there is nothing to
-        ! score and we can move on to the next tally.
-        if (filter_matches(i_filt) % bins_size() == 0) cycle TALLY_LOOP
-
-        ! Set the index of the bin used in the first filter combination
-        call filter_matches(i_filt) % set_i_bin(1)
-      end do
-
-      ! ========================================================================
-      ! Loop until we've covered all valid bins on each of the filters.
-
-      FILTER_LOOP: do
-
-        ! Reset scoring index and weight
-        filter_index = 1
-        filter_weight = ONE
-
-        ! Determine scoring index and weight for this filter combination
-        do j = 1, t % n_filters()
-          i_filt = t % filter(j)
-          i_bin = filter_matches(i_filt) % i_bin()
-          filter_index = filter_index + (filter_matches(i_filt) &
-               % bins_data(i_bin) - 1) * t % stride(j)
-          filter_weight = filter_weight * filter_matches(i_filt) &
-               % weights_data(i_bin)
-        end do
-
-        ! ======================================================================
-        ! Nuclide logic
-
-        ! Check for nuclide bins
-        NUCLIDE_LOOP: do k = 1, t % n_nuclide_bins()
-          ! Get index of nuclide in nuclides array
-          i_nuclide = t % nuclide_bins(k)
-
-          if (i_nuclide > 0) then
-            if (p % material /= MATERIAL_VOID) then
-              ! Get pointer to current material
-              mat => materials(p % material)
-
-              ! Determine index of nuclide in Material % atom_density array
-              j = mat % mat_nuclide_index(i_nuclide)
-              if (j == 0) cycle NUCLIDE_LOOP
-
-              ! Copy corresponding atom density
-              atom_density = mat % atom_density(j)
-            end if
-          else
-            atom_density = ZERO
-          end if
-
-          call score_general(p, i_tally, (k-1)*t % n_score_bins, filter_index, &
-               i_nuclide, atom_density, filter_weight)
-        end do NUCLIDE_LOOP
-
-        ! ======================================================================
-        ! Filter logic
-
-        ! Increment the filter bins, starting with the last filter to find the
-        ! next valid bin combination
-        finished = .true.
-        do j = t % n_filters(), 1, -1
-          i_filt = t % filter(j)
-          if (filter_matches(i_filt) % i_bin() < filter_matches(i_filt) % &
-               bins_size()) then
-            call filter_matches(i_filt) % set_i_bin(filter_matches(i_filt) % i_bin() + 1)
-            finished = .false.
-            exit
-          else
-            call filter_matches(i_filt) % set_i_bin(1)
-          end if
-        end do
-
-        ! Once we have finished all valid bins for each of the filters, exit
-        ! the loop.
-        if (finished) exit FILTER_LOOP
-
-      end do FILTER_LOOP
-
-      ! If the user has specified that we can assume all tallies are spatially
-      ! separate, this implies that once a tally has been scored to, we needn't
-      ! check the others. This cuts down on overhead when there are many
-      ! tallies specified
-
-      if (assume_separate) exit TALLY_LOOP
-
-      end associate
-    end do TALLY_LOOP
-
-    ! Reset filter matches flag
-    filter_matches(:) % bins_present = .false.
-
-  end subroutine score_analog_tally_mg
-
 !===============================================================================
 ! SCORE_FISSION_EOUT handles a special case where we need to store neutron
 ! production rate with an outgoing energy filter (think of a fission matrix). In
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index 97793a070d..104885edd7 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -143,7 +143,7 @@ adaptor_type<3> tally_results(int idx)
   return xt::adapt(results, size, xt::no_ownership(), shape);
 }
 
-//! Score tallies based on a simple count of events.
+//! Score tallies based on a simple count of events (for continuous energy).
 //
 //! Analog tallies ar etriggered at every collision, not every event.
 
@@ -210,7 +210,9 @@ score_analog_tally_ce(Particle* p)
           auto i_nuclide = tally.nuclides_[i];
 
           // Tally this event in the present nuclide bin if that bin represents
-          // the event nuclide or the total material.
+          // the event nuclide or the total material.  Note that the i_nuclide
+          // and flux arguments for score_general are not used for analog
+          // tallies.
           if (i_nuclide == p->event_nuclide || i_nuclide == -1)
             score_general_ce(p, i_tally, i*n_score_bins, filter_index,
               -1, -1., filter_weight);
@@ -219,8 +221,6 @@ score_analog_tally_ce(Particle* p)
         // In the case that the user has requested to tally all nuclides, we
         // can take advantage of the fact that we know exactly how nuclide
         // bins correspond to nuclide indices.  First, tally the nuclide.
-        // Note that the i_nuclide and flux arguments for score_general are
-        // not used for analog tallies.
         auto i = p->event_nuclide;
         score_general_ce(p, i_tally, i*n_score_bins, filter_index,
           -1, -1., filter_weight);
@@ -266,6 +266,117 @@ next_tally:
     match.bins_present_ = false;
 }
 
+//! Score tallies based on a simple count of events (for multigroup).
+//
+//! Analog tallies ar etriggered at every collision, not every event.
+
+extern "C" void
+score_analog_tally_mg(Particle* p)
+{
+  for (auto i_tally : model::active_analog_tallies) {
+    //TODO: off-by-one
+    const Tally& tally {*model::tallies[i_tally-1]};
+    auto n_score_bins = tally_get_n_score_bins(i_tally);
+
+    //--------------------------------------------------------------------------
+    // Loop through all relevant filters and find the filter bins and weights
+    // for this event.
+
+    // Find all valid bins in each filter if they have not already been found
+    // for a previous tally.
+    for (auto i_filt : tally.filters()) {
+      //TODO: off-by-one
+      auto& match {simulation::filter_matches[i_filt-1]};
+      if (!match.bins_present_) {
+        match.bins_.clear();
+        match.weights_.clear();
+        //TODO: off-by-one
+        model::tally_filters[i_filt-1]
+          ->get_all_bins(p, tally.estimator_, match);
+        match.bins_present_ = true;
+      }
+
+      // If there are no valid bins for this filter, then there is nothing to
+      // score so we can move on to the next tally.
+      if (match.bins_.size() == 0) goto next_tally;
+
+      // Set the index of the bin used in the first filter combination
+      match.i_bin_ = 1;
+    }
+
+    //--------------------------------------------------------------------------
+    // Loop over filter bins and nuclide bins
+
+    for (bool filter_loop_done = false; !filter_loop_done; ) {
+
+      //------------------------------------------------------------------------
+      // Filter logic
+
+      // Determine scoring index and weight for the current filter combination
+      int filter_index = 1;
+      double filter_weight = 1.;
+      for (auto i = 0; i < tally.filters().size(); ++i) {
+        auto i_filt = tally.filters(i);
+        //TODO: off-by-one
+        auto& match {simulation::filter_matches[i_filt-1]};
+        auto i_bin = match.i_bin_;
+        //TODO: off-by-one
+        filter_index += (match.bins_[i_bin-1] - 1) * tally.strides(i);
+        filter_weight *= match.weights_[i_bin-1];
+      }
+
+      //------------------------------------------------------------------------
+      // Nuclide logic
+
+      for (auto i = 0; i < tally.nuclides_.size(); ++i) {
+        auto i_nuclide = tally.nuclides_[i];
+
+        double atom_density = 0.;
+        if (i_nuclide > 0) {
+          auto j = material_nuclide_index(p->material, i_nuclide);
+          if (j == 0) continue;
+          atom_density = material_atom_density(p->material, j);
+        }
+
+        score_general_mg(p, i_tally, i*n_score_bins, filter_index,
+          i_nuclide, atom_density, filter_weight);
+      }
+
+      //------------------------------------------------------------------------
+      // Further filter logic
+
+      // Increment the filter bins, starting with the last filter to find the
+      // next valid bin combination
+      filter_loop_done = true;
+      for (int i = tally.filters().size()-1; i >= 0; --i) {
+        auto i_filt = tally.filters(i);
+        //TODO: off-by-one
+        auto& match {simulation::filter_matches[i_filt-1]};
+        if (match.i_bin_ < match.bins_.size()) {
+          ++match.i_bin_;
+          filter_loop_done = false;
+          break;
+        } else {
+          match.i_bin_ = 1;
+        }
+      }
+    }
+
+    // If the user has specified that we can assume all tallies are spatially
+    // separate, this implies that once a tally has been scored to, we needn't
+    // check the others. This cuts down on overhead when there are many
+    // tallies specified
+    if (settings::assume_separate) break;
+
+next_tally:
+    ;
+  }
+
+  // Reset all the filter matches for the next tally event.
+  for (auto& match : simulation::filter_matches)
+    match.bins_present_ = false;
+}
+
 //! Score tallies using a tracklength estimate of the flux.
 //
 //! This is triggered at every event (surface crossing, lattice crossing, or