diff --git a/src/EFTfitterDensity.jl b/src/EFTfitterDensity.jl
index b70c412..08261f4 100755
--- a/src/EFTfitterDensity.jl
+++ b/src/EFTfitterDensity.jl
@@ -59,10 +59,18 @@ struct HasLimits <: LimitsIndicator end
 struct NoLimits <: LimitsIndicator end
 # TODO: Limits not yet used
 
+#------------- Bounds Check Indicator ---------------------------------------------------------#
+# Status indicating whether observable bounds should be checked or not
+abstract type BoundsCheckIndicator end
+struct BoundsCheck <: BoundsCheckIndicator end
+struct NoBoundsCheck <: BoundsCheckIndicator end
+
+
 
 #------------- EFTfitterDensity ---------------------------------------------------------#
 struct EFTfitterDensity{
-    M<:CovOrInverseMatrix, 
+    M<:CovOrInverseMatrix,
+    B<:BoundsCheckIndicator, 
     MU<:ModelUncertaintiesIndicator, 
     NC<:NuisanceCorrelationsIndicator, 
     L<:LimitsIndicator
@@ -73,7 +81,7 @@ struct EFTfitterDensity{
     observable_maxs::Vector{Float64}
     observable_weights::Vector{Float64} # observable weights 
     crossmatrix::M # CovMatrix or InvCovMatrix
-    check_bounds::Bool  # whether to check observable bounds or not
+    check_bounds::B # whether to check observable bounds or not
     predictions::Matrix{Float64} 
     prediction_uncertainties::Matrix{Float64} # only used if model uncertainties are present
     # limit_distributions::Vector{Distribution} # only used if limits are present
@@ -128,7 +136,7 @@ function EFTfitterDensity(m::EFTfitterModel)
 
     upper_bounds = any(x->x!=Inf, observable_maxs)
     lower_bounds = any(x->x!=-Inf, observable_mins)
-    check_bounds = any([upper_bounds, lower_bounds])
+    check_bounds = any([upper_bounds, lower_bounds]) ? BoundsCheck() : NoBoundsCheck()
 
     # check if model uncertainties are present
     #TODO: make this a function
diff --git a/src/EFTfitterLikelihood.jl b/src/EFTfitterLikelihood.jl
index 0e6a516..8a00380 100644
--- a/src/EFTfitterLikelihood.jl
+++ b/src/EFTfitterLikelihood.jl
@@ -3,9 +3,14 @@ function iswithinbounds(r::Float64, min::Float64, max::Float64)
 end
 
 # TODO: add dispatch, return 1 by default
-function check_obs_bounds(r::Vector{Float64}, mins::Vector{Float64}, maxs::Vector{Float64})
-    withinbounds = [iswithinbounds(r[i], mins[i], maxs[i]) for i in 1:length(r)]
-    all(withinbounds) ? (return 1.) : (return -Inf)
+function check_observable_bounds(cb::NoBoundsCheck, predictions, mins::Vector{Float64}, maxs::Vector{Float64})
+    return 1.0
+end
+
+function check_observable_bounds(cb::BoundsCheck, predictions, mins::Vector{Float64}, maxs::Vector{Float64})
+    r = view(predictions, Threads.threadid(), :)
+    withinbounds = [iswithinbounds(r[i], mins[i], maxs[i]) for i in eachindex(r)]
+    all(withinbounds) ? (return 1.) : (return 1e50)
 end
 
 
@@ -18,7 +23,7 @@ function evaluate_funcs!(nomuncs::NoModelUncertainties, D::EFTfitterDensity, cur
     end
 end
 
-# with model uncertainties #TODO: remove arr 
+# with model uncertainties
 function evaluate_funcs!(muncs::HasModelUncertainties, D::EFTfitterDensity, current_params)
     funcs = D.observable_functions
     for i in eachindex(funcs)
@@ -35,10 +40,11 @@ function DensityInterface.logdensityof(
 )
     evaluate_funcs!(D.model_uncertainties, D, params)
 
-    #TODO: check_observable_bounds
+    bounds_factor = check_observable_bounds(D.check_bounds, D.predictions, D.observable_mins, D.observable_maxs)
+
     result = calculate_likelihood(D.model_uncertainties, D.nuisance_correlations, D, params)
 
-    return result
+    return bounds_factor * result
 end
 
 # no model uncertainties & no limits, weights are already included in the inverse covariance matrix
diff --git a/test/test_plain.jl b/test/test_plain.jl
index 3ff2a03..fba398e 100644
--- a/test/test_plain.jl
+++ b/test/test_plain.jl
@@ -23,7 +23,7 @@ measurements = (
     meas2 = Measurement(Observable(testfunc1, min=0, max=1000), 222.2, 
     uncertainties = (unc1=20.1, unc2=20.2, unc3=23.3), active=false),
     
-    meas3 = Measurement(Observable(testfunc1, min=0, max=1000), 333.3, 
+    meas3 = Measurement(Observable(testfunc1, min=0, max=25000), 333.3, 
     uncertainties = (unc1=30.1, unc2=30.2, unc3=30.3), active=true),
     
     meas4 = MeasurementDistribution(Function[testfunc1, testfunc1, testfunc1],
@@ -54,12 +54,12 @@ correlations = (
 )
 
 
-model = EFTfitterModel(parameters, measurements, correlations)
+model = EFTfitterModel(parameters, measurements, correlations, )
 posterior = PosteriorMeasure(model)
 
 v = (p1 = 10.826122384321511, p2 = -8.32129957354641)
 logp = logdensityof(posterior)
-
+logp(v)
 @test logp(v) ≈ -5.2063526518e9
 
 
@@ -72,6 +72,6 @@ t = @benchmark logp(v)
 
 # with @SVector
 t = @benchmark logp(v)
-@test t.allocs == 11
-@test t.memory == 400
-@test minimum(t.times) ≈ 163 atol=5
+@test t.allocs == 12
+@test t.memory == 464
+@test minimum(t.times) ≈ 191 atol=5