diff --git a/resolve/likelihood.py b/resolve/likelihood.py
index 99e326433060b27e756ab5c333102394090c6bfc..885f2198f3e168306b6fafb45450dcb999ee591a 100644
--- a/resolve/likelihood.py
+++ b/resolve/likelihood.py
@@ -198,11 +198,10 @@ def CalibrationLikelihood(
     cops = _duplicate(_obj2list(calibration_operator, ift.Operator), len(obs))
     log_icovs = _duplicate(_obj2list(log_inverse_covariance_operator, ift.Operator), len(obs))
     model_d = _duplicate(_obj2list(model_visibilities, ift.Field), len(obs))
-    model_d = [ift.makeOp(mm) @ cop for mm, cop in zip(model_d, cops)]
 
     if len(obs) > 1:
         raise NotImplementedError
-    obs, model_d, log_icov = obs[0], model_d[0], log_icovs[0]
+    obs, model_d, log_icov, cop = obs[0], model_d[0], log_icovs[0]; cops[0]
 
     dt = DtypeConverter(model_d.target, np.complex128, obs.vis.dtype)
     dt_icov = DtypeConverter(model_d.target, np.float64, obs.weight.dtype)
@@ -213,11 +212,12 @@ def CalibrationLikelihood(
             inverse_covariance=obs.weight,
             nthreads=nthreads,
             mask=obs.mask,
+            multiplicative=dt(model_d)
         )
-        return e @ dt @ model_d
+        return e @ dt @ cop
     else:
         s0, s1 = "model data", "inverse covariance"
         e = VariableCovarianceDiagonalGaussianLikelihood(
             obs.vis, s0, s1, mask=obs.mask, nthreads=nthreads
         )
-        return e @ ((dt @ model_d).ducktape_left(s0) + (dt_icov @ log_icov).ducktape_left(s1))
+        return e @ ((dt @ model_d @ cop).ducktape_left(s0) + (dt_icov @ log_icov).ducktape_left(s1))