diff --git a/cli/dumpTree.py b/cli/dumpTree.py
index b1527214..280f5adc 100644
--- a/cli/dumpTree.py
+++ b/cli/dumpTree.py
@@ -23,6 +23,7 @@
                            ("dx", "f4"), ("long_diff", "f4"),
                            ("pixel_plane", "i4"), ("t_end", "f4"),
                            ("dEdx", "f4"), ("dE", "f4"), ("t", "f4"),
+                           ("dEdx_nonIonizing", "f4"), ("dE_nonIonizing", "f4"),
                            ("y", "f4"), ("x", "f4"), ("z", "f4"),
                            ("n_photons","f4")], align=True)
 
@@ -302,6 +303,7 @@ def dump(input_file, output_file):
                 segment[iHit]["t0_end"] = hitSegment.GetStop().T() * edep2us
                 segment[iHit]["t_end"] = 0
                 segment[iHit]["dE"] = hitSegment.GetEnergyDeposit()
+                segment[iHit]["dE_nonIonizing"] = hitSegment.GetSecondaryDeposit()
                 xd = segment[iHit]["x_end"] - segment[iHit]["x_start"]
                 yd = segment[iHit]["y_end"] - segment[iHit]["y_start"]
                 zd = segment[iHit]["z_end"] - segment[iHit]["z_start"]
@@ -313,6 +315,7 @@ def dump(input_file, output_file):
                 segment[iHit]["t0"] = (segment[iHit]["t0_start"] + segment[iHit]["t0_end"]) / 2.
                 segment[iHit]["t"] = 0
                 segment[iHit]["dEdx"] = hitSegment.GetEnergyDeposit() / dx if dx > 0 else 0
+                segment[iHit]["dEdx_nonIonizing"] = hitSegment.GetSecondaryDeposit() / dx if dx > 0 else 0
                 segment[iHit]["pdgId"] = trajectories[hitSegment.Contrib[0]]["pdgId"]
                 segment[iHit]["n_electrons"] = 0
                 segment[iHit]["long_diff"] = 0
diff --git a/larndsim/quenching.py b/larndsim/quenching.py
index 6081e3ec..d3f934ef 100644
--- a/larndsim/quenching.py
+++ b/larndsim/quenching.py
@@ -26,19 +26,25 @@ def quench(tracks, mode):
         dEdx = tracks[itrk]["dEdx"]
         dE = tracks[itrk]["dE"]
 
+        dEdx_nonIonizing = tracks[itrk]["dEdx_secondary"]
+        dE_nonIonizing = tracks[itrk]["dE_secondary"]
+
+        dEdx_ionizing = dEdx - dEdx_nonIonizing
+        dE_ionizing = dE - dE_nonIonizing
+
         recomb = 0
         if mode == physics.BOX:
             # Baller, 2013 JINST 8 P08005
-            csi = physics.BOX_BETA * dEdx / (detector.E_FIELD * detector.LAR_DENSITY)
+            csi = physics.BOX_BETA * dEdx_ionizing / (detector.E_FIELD * detector.LAR_DENSITY)
             recomb = max(0, log(physics.BOX_ALPHA + csi)/csi)
         elif mode == physics.BIRKS:
             # Amoruso, et al NIM A 523 (2004) 275
-            recomb = physics.BIRKS_Ab / (1 + physics.BIRKS_kb * dEdx / (detector.E_FIELD * detector.LAR_DENSITY))
+            recomb = physics.BIRKS_Ab / (1 + physics.BIRKS_kb * dEdx_ionizing / (detector.E_FIELD * detector.LAR_DENSITY))
         else:
             raise ValueError("Invalid recombination mode: must be 'physics.BOX' or 'physics.BIRKS'")
 
         if isnan(recomb):
             raise RuntimeError("Invalid recombination value")
 
-        tracks[itrk]["n_electrons"] = recomb * dE / physics.W_ION
+        tracks[itrk]["n_electrons"] = recomb * dE_ionizing / physics.W_ION
         tracks[itrk]["n_photons"] = (dE/light.W_PH - tracks[itrk]["n_electrons"]) * light.SCINT_PRESCALE