diff --git a/src/constants.F90 b/src/constants.F90
index a21a9c5915..4186ebfb11 100644
--- a/src/constants.F90
+++ b/src/constants.F90
@@ -252,7 +252,7 @@ module constants
        EVENT_FISSION =  3 
 
   ! Tally score type
-  integer, parameter :: N_SCORE_TYPES = 17
+  integer, parameter :: N_SCORE_TYPES = 19
   integer, parameter :: &
        SCORE_FLUX       = -1,  & ! flux
        SCORE_TOTAL      = -2,  & ! total reaction rate
@@ -269,8 +269,10 @@ module constants
        SCORE_ABSORPTION = -13, & ! absorption rate
        SCORE_FISSION    = -14, & ! fission rate
        SCORE_NU_FISSION = -15, & ! neutron production rate
-       SCORE_CURRENT    = -16, & ! partial current
-       SCORE_EVENTS     = -17    ! number of events
+       SCORE_K_FISSION  = -16, & ! fission energy production rate
+       SCORE_CHI        = -17, & ! outgoing fission neutron energy distribution
+       SCORE_CURRENT    = -18, & ! partial current
+       SCORE_EVENTS     = -19    ! number of events
        
   ! Maximum scattering order supported
   integer, parameter :: SCATT_ORDER_MAX = 10
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 4bfa657855..e67073508d 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -1853,6 +1853,16 @@ contains
             end if
           case ('nu-fission')
             t % score_bins(j) = SCORE_NU_FISSION
+          case ('k-fission')
+            t % score_bins(j) = SCORE_K_FISSION
+            
+            ! Set tally estimator to analog
+            t % estimator = ESTIMATOR_ANALOG
+          case ('chi')
+            t % score_bins(j) = SCORE_CHI
+            
+            ! Set tally estimator to analog
+            t % estimator = ESTIMATOR_ANALOG
           case ('current')
             t % score_bins(j) = SCORE_CURRENT
             t % type = TALLY_SURFACE_CURRENT
diff --git a/src/output.F90 b/src/output.F90
index 8af28c1b59..1dae48a5cc 100644
--- a/src/output.F90
+++ b/src/output.F90
@@ -821,6 +821,10 @@ contains
         string = trim(string) // ' fission'
       case (SCORE_NU_FISSION)
         string = trim(string) // ' nu-fission'
+      case (SCORE_K_FISSION)
+        string = trim(string) // ' k-fission'
+      case (SCORE_CHI)
+        string = trim(string) // ' chi'
       case (SCORE_CURRENT)
         string = trim(string) // ' current'
       end select
@@ -1461,6 +1465,8 @@ contains
     score_names(abs(SCORE_ABSORPTION)) = "Absorption Rate"
     score_names(abs(SCORE_FISSION))    = "Fission Rate"
     score_names(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate"
+    score_names(abs(SCORE_K_FISSION)) = "Kappa-Fission Rate"
+    score_names(abs(SCORE_CHI)) = "Chi"
     score_names(abs(SCORE_EVENTS))     = "Events"
 
     ! Create filename for tally output
diff --git a/src/tally.F90 b/src/tally.F90
index 8e90010db9..ffafcdbd4b 100644
--- a/src/tally.F90
+++ b/src/tally.F90
@@ -315,6 +315,26 @@ contains
               score = keff * p % wgt_bank
 
             end if
+          
+          case (SCORE_K_FISSION)
+            ! Skip any non-fission events
+            if (p % event /= EVENT_FISSION) cycle SCORE_LOOP
+
+            ! All fission events will contribute, so again we can use
+            ! particle's weight entering the collision as the estimate for
+            ! the fission energy production rate
+            ! The ENDF standard (ENDF-102) states that MT 18 stores
+            ! the fission energy as the Q_value, so this will be obtained from
+            ! the MT of index_fission(1); the standard is mute on
+            ! whether or not the other fission channels store their own
+            ! fission energy output; after inspection of the ACE files, all
+            ! nuclides observed with multiple fission channels have the same 
+            ! Q-value for all of them; therefore, this code just uses the first
+            ! fission reaction instead of storing the actual fission type which
+            ! occured.
+            n = nuclides(p % event_nuclide) % index_fission(1)
+            score = last_wgt * &
+              nuclides(p % event_nuclide) % reactions(n) % Q_value
 
           case (SCORE_EVENTS)
             ! Simply count number of scoring events