diff --git a/R/baseline-GCIMSDataset.R b/R/baseline-GCIMSDataset.R
index 8910384..cc5d04c 100644
--- a/R/baseline-GCIMSDataset.R
+++ b/R/baseline-GCIMSDataset.R
@@ -11,12 +11,12 @@
 setMethod(
   "estimateBaseline",
   "GCIMSDataset",
-  function(object, dt_peak_fwhm_ms, dt_region_multiplier, rt_length_s) {
+  function(object, dt_peak_fwhm_ms, dt_region_multiplier, rt_length_s, remove =TRUE) {
 
     delayed_op <- DelayedOperation(
       name = "estimateBaseline",
       fun = estimateBaseline,
-      params = list(dt_peak_fwhm_ms = dt_peak_fwhm_ms, dt_region_multiplier = dt_region_multiplier, rt_length_s = rt_length_s)
+      params = list(dt_peak_fwhm_ms = dt_peak_fwhm_ms, dt_region_multiplier = dt_region_multiplier, rt_length_s = rt_length_s, remove = remove)
     )
     object$appendDelayedOp(delayed_op)
     # We recompute these, but  maybe we could just reset them to zero...
diff --git a/R/baseline-GCIMSSample.R b/R/baseline-GCIMSSample.R
index f818f07..0d2d6b3 100644
--- a/R/baseline-GCIMSSample.R
+++ b/R/baseline-GCIMSSample.R
@@ -15,7 +15,7 @@
 #' @export
 methods::setMethod(
   "estimateBaseline", "GCIMSSample",
-  function(object, dt_peak_fwhm_ms, dt_region_multiplier, rt_length_s) {
+  function(object, dt_peak_fwhm_ms, dt_region_multiplier, rt_length_s, remove) {
     rt <- rtime(object)
     dt <- dtime(object)
     int <- intensity(object)
@@ -28,6 +28,9 @@ methods::setMethod(
     rt_basel <- estimate_baseline_tr(int - dt_basel, region_size = rt_region_size_pts)
     full_basel <- rt_basel + dt_basel
     baseline(object) <- full_basel
+    if (remove){
+      object@data <- int - full_basel
+    }
     object
   }
 )