diff --git a/src/treeMeasure.jl b/src/treeMeasure.jl
index ec99880..aee82c9 100644
--- a/src/treeMeasure.jl
+++ b/src/treeMeasure.jl
@@ -121,15 +121,24 @@ end
 """
     measure2siteoperator(net::TreeNetwork, O1, O2, site1::Int64, site2::Int64})
 
-For a Tree, compute the local expectation value of two one-site operators O1 and O2 only for site1 and site2.
+For a Tree, compute the local expectation value of two one-site operators O1 and O2 only for site1 and site2. If O1 can only be applied on site1 and O2 on site2, compute the local expectation value of the one-site operator O1 on site1 and of the one-site operator O2 on site2
 
 """
+
 function measure2siteoperator(net::TreeNetwork, O1, O2, site1::Int64,site2::Int64)
 
     sites=(site1,site2)
-    herm_cis = ishermitian(O1*O2)
-    herm_trans = ishermitian(O1) && ishermitian(O2)
-    conjpair = O1 == O2'
+    if size(O1)==size(O2)
+    	herm_cis = ishermitian(O1*O2)
+    	herm_trans = ishermitian(O1) && ishermitian(O2)
+    	conjpair = O1 == O2'
+	T = (herm_cis && herm_trans) ? Float64 : ComplexF64
+        expval = Array{T,2}(undef, 2,2)
+    else
+	herm_trans = ishermitian(O1) && ishermitian(O2)
+        T = herm_trans ? Float64 : ComplexF64
+        expval = Vector{T}(undef, 1)
+    end
 
     if site1==site2
         return measure1siteoperator(A, O1*O2, site1)
@@ -141,44 +150,57 @@ function measure2siteoperator(net::TreeNetwork, O1, O2, site1::Int64,site2::Int6
     firstsite = sites[1]
     lastsite = sites[2]
     N = length(path(net, firstsite, lastsite))
-    T = (herm_cis && herm_trans) ? Float64 : ComplexF64
-    expval = Array{T,2}(undef, 2,2)
 
     if !in(findheadnode(net), path(net, firstsite, lastsite)[2:end-1])
         if length(path(net, firstsite)) < length(path(net, lastsite))
             nearestsite = firstsite
             farthestsite = lastsite
+	    m1 = O1
+	    m2 = O2
         else
             nearestsite = lastsite
             farthestsite = firstsite
+            m1 = O2
+            m2 = O1
         end
         for (A, id, dir) in Path(net, nearestsite)[1:end-1]
             ρ = rhoAAstar(ρ, A, 1, dir)
         end
         for (i, (A, id, dir)) in enumerate(Path(net, nearestsite, farthestsite)[1:end-1])
-	    if id == nearestsite
+	    if id == nearestsite && size(O1)==size(O2)
                v = rhoAOAstar(ρ, A, O1*O2, 1, nothing)
                expval[1,1] = v
 	    end
-            ρ1 = rhoAOAstar(ρ, A, O1, 1, dir)
-            ρ2 = rhoAOAstar(ρ, A, O2, 1, dir)
+	    if size(O1)==size(O2)
+               ρ1 = rhoAOAstar(ρ, A, O1, 1, dir)
+               ρ2 = rhoAOAstar(ρ, A, O2, 1, dir)
+	    else
+	       ρ1 = rhoAOAstar(ρ, A, m1, 1, dir)
+	    end
             for (j, (B, jd, djr)) in enumerate(Path(net, id, farthestsite)[2:end])
-		if id==nearestsite
+		if id==nearestsite && size(O1)==size(O2)
                    v = rhoAOAstar(ρ1, B, O2, 1, nothing)
                    expval[1,2] = v
                    v = rhoAOAstar(ρ2, B, O1, 1, nothing)
                    expval[2,1] = v
+		elseif id==nearestsite && size(O1)!=size(O2)
+		   expval = rhoAOAstar(ρ1, B, m2, 1, nothing)
+		   T<:Real && (expval=real(expval))
 		end
                 djr != nothing && (ρ1 = rhoAAstar(ρ1, B, 1, djr))
-                djr != nothing && (ρ2 = rhoAAstar(ρ2, B, 1, djr))
+		if size(O1)==size(O2)
+                   djr != nothing && (ρ2 = rhoAAstar(ρ2, B, 1, djr))
+		end
             end
             ρ = rhoAAstar(ρ, A, 1, dir)
         end
-        v = rhoAOAstar(ρ, net[farthestsite], O1*O2, 1, nothing)
-        expval[2,2] = v
-        if nearestsite == lastsite
-            expval = reverse(reverse(expval, dims=1),dims=2)
-        end
+	if size(O1)==size(O2)
+           v = rhoAOAstar(ρ, net[farthestsite], O1*O2, 1, nothing)
+           expval[2,2] = v
+           if nearestsite == lastsite
+               expval = reverse(reverse(expval, dims=1),dims=2)
+           end
+	end
         return expval
     else
         throw("missing method")