diff --git a/src/LorentzVectorHEP.jl b/src/LorentzVectorHEP.jl
index 4504f76..1159930 100644
--- a/src/LorentzVectorHEP.jl
+++ b/src/LorentzVectorHEP.jl
@@ -4,14 +4,22 @@ using LorentzVectors # provides x, y, z, t
 
 export LorentzVectorCyl, LorentzVector
 
-export px, py, pz, energy, fast_mass, pt, rapidity, eta, phi, mass
-export deltaphi, deltar, deltaeta
+
+export px, py, pz, energy, fast_mass, pt, rapidity, eta, phi, theta, mass
+export deltaphi, deltar, deltaeta, deltatheta
 export ΔR, Δϕ, Δη
 export fromPtEtaPhiE
 
 include("cartesian.jl")
 include("cylindrical.jl")
 
+const Δϕ = deltaphi
+const Δη = deltaeta
+const ΔR = deltar
+const Δθ = deltatheta
+
+export ΔR, Δϕ, Δη, Δθ
+
 # conversion
 function LorentzVector(v::LorentzVectorCyl)
     x = px(v)
diff --git a/src/cartesian.jl b/src/cartesian.jl
index 720d0df..94629ca 100644
--- a/src/cartesian.jl
+++ b/src/cartesian.jl
@@ -52,6 +52,14 @@ function phi02pi(lv::LorentzVector)
     return phi(lv) < 0.0 ? phi(lv) + 2π : phi(lv)
 end
 
+function theta(lv::LorentzVector)
+    return lv.x == 0.0 && lv.y == 0.0 && lv.z == 0.0 ? 0.0 : atan2(pt(lv), lv.z)
+end
+
+function deltatheta(lv1::LorentzVector, lv2::LorentzVector)
+    return phi_mpi_pi(theta(lv1) - theta(lv2))
+end
+
 function phi_mpi_pi(x)
     twopi = 2pi
     while (x >= pi)
diff --git a/src/cylindrical.jl b/src/cylindrical.jl
index b1ead7c..9dcbab9 100644
--- a/src/cylindrical.jl
+++ b/src/cylindrical.jl
@@ -89,12 +89,10 @@ end
     dphi = dphi - 2*pi*(dphi > pi) + 2*pi*(dphi <= -pi)
     return dphi
 end
-const Δϕ = deltaphi
 
 @inline function deltaeta(v1::LorentzVectorCyl, v2::LorentzVectorCyl)
     return v2.eta - v1.eta
 end
-const Δη = deltaeta
 
 @inline function deltar2(v1::LorentzVectorCyl, v2::LorentzVectorCyl)
     dphi = deltaphi(v1,v2)
@@ -102,12 +100,11 @@ const Δη = deltaeta
     return muladd(dphi, dphi, deta^2)
 end
 deltar(v1::LorentzVectorCyl, v2::LorentzVectorCyl) = sqrt(deltar2(v1, v2))
-const ΔR = deltar
 
+const ΔR = deltar
 
 function fromPtEtaPhiE(pt, eta, phi, E) 
     m2 = E^2 - pt^2 - (sinh(eta) * pt)^2
     m = sign(m2) * sqrt(abs(m2)) 
     return LorentzVectorCyl(pt, eta, phi, m) 
-end
-  
+end
\ No newline at end of file