diff --git a/4_grains/gtComputeIncomingBeamIntensity.m b/4_grains/gtComputeIncomingBeamIntensity.m
index 47080eec83e1c57f909939587d88a143e0070dd5..8e79e581768be4e15af9002c99d7d47bec6bac1d 100644
--- a/4_grains/gtComputeIncomingBeamIntensity.m
+++ b/4_grains/gtComputeIncomingBeamIntensity.m
@@ -7,7 +7,14 @@ function [gr, refs] = gtComputeIncomingBeamIntensity(gr, parameters, refs)
     end
     omega_step = 180 / parameters.acq.nproj;
 
-    gc_det_pos = gtGeoGrainCenterSam2Det(gr.center, gr.allblobs.omega, parameters);
+    if (isfield(gr.proj, 'ondet'))
+        gr_included = gr.proj.ondet(gr.proj.included);
+    else
+        gr_included = gr.ondet(gr.included);
+    end
+
+    omegas = gr.allblobs.omega(gr_included);
+    gc_det_pos = gtGeoGrainCenterSam2Det(gr.center, omegas, parameters);
     rounded_gc_det = round(gc_det_pos);
 
     img_nums = omegas / omega_step;
@@ -24,14 +31,14 @@ function [gr, refs] = gtComputeIncomingBeamIntensity(gr, parameters, refs)
 
     max_selection = max_cs > 1e-4;
     min_indx = sub2ind(size(refs), rounded_gc_det(:, 2), rounded_gc_det(:, 1), min_ws);
-    max_indx = sub2ind(size(refs), rounded_gc_det(:, 2), rounded_gc_det(:, 1), max_ws(max_selection));
+    max_indx = sub2ind(size(refs), rounded_gc_det(max_selection, 2), rounded_gc_det(max_selection, 1), max_ws(max_selection));
     min_vals = refs(min_indx) .* min_cs;
     max_vals = refs(max_indx) .* max_cs(max_selection);
 
     int_profile = min_vals;
     int_profile(max_selection) = int_profile(max_selection) + max_vals;
 
-    gr.allblobs.beam_intensity = int_profile;
+    gr.beam_intensity = int_profile;
 
 %     ref0 = refs;
 %     ref0(min_indx) = -20000;