diff --git a/src/BPGates.jl b/src/BPGates.jl
index ae786c5..8888d70 100644
--- a/src/BPGates.jl
+++ b/src/BPGates.jl
@@ -558,41 +558,41 @@ function QuantumClifford.apply!(state::BellState, g::T1NoiseOp)
     output_state = if input_state==1
         if     r < 0.5*λ₁^2 - λ₁ + 1
             1
-        elseif r < 0.5*λ₁^2 - λ₁ + 1  +  0.5*λ₁*(1-λ₁)
-            2
-        elseif r < 0.5*λ₁^2 - λ₁ + 1  +  0.5*λ₁*(1-λ₁)  +  0.5*λ₁^2
-            3
+        elseif r < 0.5*λ₁^2 - λ₁ + 1  +  0.5*λ₁^2 
+            2 # XXX
+        elseif r < 0.5*λ₁^2 - λ₁ + 1  +  0.5*λ₁^2 +  0.5*λ₁*(1-λ₁)  
+            3 # XXX
         else # r < 1 = 0.5*λ₁^2 - λ₁ + 1  +  0.5*λ₁*(1-λ₁)  +  0.5*λ₁^2  +   0.5*λ₁*(1-λ₁)
             4
         end
     elseif input_state==2
-        if     r < 0.5*λ₁
+        if     r < 0.5*λ₁^2
             1
-        elseif r < 0.5*λ₁ + 1-λ₁
-            2
-        elseif r < 0.5*λ₁ + 1-λ₁  +  0.5*λ₁
-            3
-        else # r < 1 = 0.5*λ₁ + 1-λ₁  +  0.5*λ₁  + 0
+        elseif r < 0.5*λ₁^2 + 0.5*λ₁^2 -λ₁ + 1
+            2 # XXX
+        elseif r < 0.5*λ₁^2 + 0.5*λ₁^2 -λ₁ + 1 + 0.5*λ₁*(1-λ₁)
+            3 # XXX
+        else # r < 1 = 0.5*λ₁^2 + 0.5*λ₁^2 -λ₁ + 1 + 0.5*λ₁*(1-λ₁) + 0.5*λ₁*(1-λ₁)
             4
         end
     elseif input_state==3
-        if     r < 0.25*λ₁*(λ₁+1)
+        if     r < 0.5*λ₁
             1
-        elseif r < 0.25*λ₁*(λ₁+1)  +  0.5*λ₁*(1-λ₁)
-            2
-        elseif r < 0.25*λ₁*(λ₁+1)  +  0.5*λ₁*(1-λ₁)  +  0.25*λ₁*(λ₁+1)
-            3
-        else # r < 1 = 0.25*λ₁*(λ₁+1)  +  0.5*λ₁*(1-λ₁)  +  0.25*λ₁*(λ₁+1)  +  0.5*λ₁*(1-λ₁)
+        elseif r < 0.5*λ₁  +  0.5*λ₁
+            2 # XXX
+        elseif r < 0.5*λ₁  +  0.5*λ₁  +  (1-λ₁)
+            3 # XXX
+        else # r < 1 = 0.5*λ₁  +  0.5*λ₁  +  (1-λ₁)  +  0
             4
         end
     else # input_state==4
         if     r < 0.5*λ₁
             1
-        # elseif r < 0.5*λ₁ + 0             # output_state 2 is never reached
-        #     2
-        elseif r < 0.5*λ₁ + 0  +  0.5*λ₁
-            3
-        else # r < 1 = 0.5*λ₁ + 0  +  0.5*λ₁ + 1-λ₁
+        elseif r < 0.5*λ₁ + 0.5*λ₁
+            2 # XXX
+        # elseif r < 0.5*λ₁ + 0.5*λ₁ + 0             # output_state 3 is never reached
+        #     3 # XXX
+        else # r < 1 = 0.5*λ₁ + 0.5*λ₁ + 0 + 1-λ₁
             4
         end
     end
diff --git a/test/test_quantumoptics.jl b/test/test_quantumoptics.jl
index fbffa22..0caf081 100644
--- a/test/test_quantumoptics.jl
+++ b/test/test_quantumoptics.jl
@@ -1,6 +1,8 @@
 @testitem "QuantumOptics comparisons" begin
 
-using BPGates, QuantumClifford, QuantumOptics
+using Test
+
+using BPGates, QuantumClifford, QuantumOpticsBase
 
 using BPGates: T1NoiseOp
 
@@ -9,30 +11,26 @@ using LinearAlgebra: diag
 # define some BP objects
 λ = 0.2
 N = 100000
-s = BellState(1)
 opBP = T1NoiseOp(1, λ)
 
-# compute using BP the density matrix after T1 noise
-ρBP = sum([dm(Ket(Stabilizer(apply!(s,op)))) for i in 1:N])/N
-
 # define some QO objects
 b = SpinBasis(1//2)
 l0 = spinup(b)
 l1 = spindown(b)
 l00 = l0⊗l0
-l01 = l0⊗l1
-l10 = l1⊗l0
+l01 = l1⊗l0 # XXX be careful with the reversed endiandness - this is ket [0 1 0 0]
+l10 = l0⊗l1 # XXX be careful with the reversed endiandness - this is ket [0 0 1 0]
 l11 = l1⊗l1
 bell00 = (l00+l11)/sqrt(2)
-bell01 = (l01+l10)/sqrt(2)
-bell10 = (l00-l11)/sqrt(2)
+bell10 = (l00-l11)/sqrt(2) # XXX                                # bellstateindex = 2
+bell01 = (l01+l10)/sqrt(2) # XXX                                # bellstateindex = 3
 bell11 = (l01-l10)/sqrt(2)
 #|`00`|`+XX +ZZ`|`∣00⟩+∣11⟩`|`∣++⟩+∣--⟩`|`∣i₊i₋⟩+∣i₋i₊⟩`|
-#|`01`|`+XX -ZZ`|`∣01⟩+∣10⟩`|`∣++⟩-∣--⟩`|`∣i₊i₊⟩-∣i₋i₋⟩`|
-#|`10`|`-XX +ZZ`|`∣00⟩-∣11⟩`|`∣+-⟩+∣-+⟩`|`∣i₊i₊⟩+∣i₋i₋⟩`|
+#|`10`|`-XX +ZZ`|`∣00⟩-∣11⟩`|`∣+-⟩+∣-+⟩`|`∣i₊i₊⟩+∣i₋i₋⟩`|       # be careful : bellstateindex = 2
+#|`01`|`+XX -ZZ`|`∣01⟩+∣10⟩`|`∣++⟩-∣--⟩`|`∣i₊i₊⟩-∣i₋i₋⟩`|       # be careful : bellstateindex = 3
 #|`11`|`-XX -ZZ`|`∣01⟩-∣10⟩`|`∣+-⟩-∣-+⟩`|`∣i₊i₋⟩-∣i₋i₊⟩`|
 
-# compute using QO the density matrix after T1 noise
+# T1 noise in the QO formalism
 krausOp1 = projector(l0) + √(1-λ) * projector(l1)
 krausOp2 = √(λ) * projector(l0, l1')
 id = identityoperator(b)
@@ -40,16 +38,29 @@ k1 = krausOp1⊗id
 k2 = krausOp2⊗id
 k3 = id⊗krausOp1
 k4 = id⊗krausOp2
-ρ0 = dm(bell00)
-ρ1  = k1*ρ0*k1' + k2*ρ0*k2'
-ρQO = k3*ρ1*k3' + k4*ρ1*k4'
 
 # switch to the Bell basis
 to_bell_basis = projector(l00,bell00')+projector(l01,bell01')+projector(l10,bell10')+projector(l11,bell11')
 
-ρbBP = to_bell_basis*ρBP*to_bell_basis'
-ρbQO = to_bell_basis*ρQO*to_bell_basis'
+for bellstateindex in 1:4
+
+    # compute using BP the density matrix after T1 noise
+    s = BellState(BPGates.int_to_bit(bellstateindex, Val(2)))
+    ρBP = sum([dm(Ket(Stabilizer(apply!(copy(s),opBP)))) for i in 1:N])/N
+
+    # compute using QO the density matrix after T1 noise
+    ψ = [bell00,bell10,bell01,bell11][bellstateindex]
 
-@test_broken isapprox(diag(ρbBP.data), diag(ρbQO.data), atol=10/sqrt(N))
+    @test abs(ψ' * Ket(Stabilizer(s))) ≈ 1
+    ρ0 = dm(ψ)
+    ρ1  = k1*ρ0*k1' + k2*ρ0*k2'
+    ρQO = k3*ρ1*k3' + k4*ρ1*k4'
+
+    ρbBP = to_bell_basis*ρBP*to_bell_basis'
+    ρbQO = to_bell_basis*ρQO*to_bell_basis'
+
+    @test isapprox(diag(ρbBP.data), diag(ρbQO.data), atol=10/sqrt(N))
+
+end
 
 end