diff --git a/src/Procdep/realint.f b/src/Procdep/realint.f
index 7ded2d2..b075384 100644
--- a/src/Procdep/realint.f
+++ b/src/Procdep/realint.f
@@ -733,6 +733,7 @@ c        call singcheck(qqb_QQb_g,qqb_QQb_gs,p)
 c       call singletop2_scale_setup(p)
 c       call singcheck(singletop2_real_heavy, singletop2_gs_heavy, p)
       elseif (kcase==kbq_tpq) then
+        bbpart = .false.
         bbfrac = 0._dp
 c        call singcheck(qqb_tbb_g,qqb_tbb_gs,p)
         if (includereal) call qqb_tbb_g(p,msq)
@@ -1222,6 +1223,7 @@ c--- for single top + b, make sure to use two different scales
       if (nproc==169 .or. nproc==164 .or. nproc==161 .or. nproc==166) then
         if ((j == 0 .and. abs(k) /= 5) .or. (abs(j) /= 5 .and. k ==0)) then
             bbfrac(0) = bbfrac(0) + xmsqjk
+            bbpart = .true.
         endif
       endif
 
@@ -1272,11 +1274,8 @@ c         endif
 
       ! summation and histogramming for all dipole contributions
       do nd=0,ndmax
-          if (kcase == ktopanom .or. kcase == kbq_tpq) then
-            if (xmsq(nd) /= 0._dp) then
-                bbfrac(nd) = bbfrac(nd)/xmsq(nd)
-            endif
-          endif
+          ! currentNd is required early, jet content tagging depends on it!
+          currentNd = nd
 
           xmsq(nd)=xmsq(nd)*flux*pswt/BrnRat
 
@@ -1284,9 +1283,11 @@ c         endif
               q(:,:) = ptilde(nd,:,:)
               ! check cuts
               if (incldip(nd)) incldip(nd)=includedipole(nd,q)
+
               if (incldip(nd) .eqv. .false.) then
                   call dotem(nvec,p,s)
                   xmsq(nd) = 0._dp
+                  bbfrac(nd) = 0._dp
               endif
           endif
 
@@ -1359,15 +1360,14 @@ c         endif
                       endif
                   endif
 
-                  currentNd = nd
                   call getptildejet(nd,pjet)
                   call dotem(nvec,pjet,s)
                   if (bbpart) then
                       ! separately bin bb part and other part
                       bbpart = .true.
-                      call nplotter_new(pjet, xmsq(nd)*wgt*bbfrac(nd))
+                      call nplotter_new(pjet, wgt*bbfrac(nd))
                       bbpart = .false.
-                      call nplotter_new(pjet, xmsq(nd)*wgt*(1-bbfrac(nd)))
+                      call nplotter_new(pjet, (xmsq(nd)-bbfrac(nd))*wgt)
                   else
                       call nplotter_new(pjet,xmsq(nd)*wgt)
                   endif
diff --git a/src/Procdep/virtint.f b/src/Procdep/virtint.f
index 0f7d906..0c16582 100644
--- a/src/Procdep/virtint.f
+++ b/src/Procdep/virtint.f
@@ -11,7 +11,7 @@
         use singletop2_m, only : singletop2_z, singletop2_tree, singletop2_virt
         use singletop2_scale_m
         use singletop_virtint, only: singletop_virtintf => virtint
-        use bbfrac_m, only : bbfrac
+        use bbfrac_m, only : bbfrac, bbpart
         use MCFMStorage
         use SCET
         use m_gencuts, only : enable_reweight_user, reweight_user
@@ -993,6 +993,7 @@ c Sanity check for tiny pt(top)/mt -- otherwise numerical Gram det. problems
         ! this has been debugged, but ocasionally there are unstable points
         ! due to the strong cancellation effects with epinv /= 0
         bbfrac = 0._dp
+        bbpart = .false.
         msqv_heavy = 0._dp
         msqv_light = 0._dp
 
@@ -1014,6 +1015,8 @@ c Sanity check for tiny pt(top)/mt -- otherwise numerical Gram det. problems
             goto 999
         endif
       elseif (kcase==kbq_tpq) then
+        bbfrac = 0._dp
+        bbpart = .false.
         call bq_tpq(p,msq)
         call bq_tpq_v(p,msqv)
         call bq_tpq_z(p,z)