relativistic optical depth

import sympy
sympy.init_printing()

b = sympy.Symbol('beta')
g = sympy.Symbol('gamma')
dx = sympy.Symbol('dx')
q = sympy.Symbol('theta')
boost = sympy.Matrix([[1,0,0,0],
                      [0,g,0,-g*b],
                      [0,0,1,0],
                      [0,-g*b,0,g]])
start = sympy.Matrix([0,-dx,0,0])
end = sympy.Matrix([sympy.sin(q)*dx/(sympy.cos(q)-b),
                    b*dx/(sympy.cos(q)-b),
                    0,
                    dx/(sympy.cos(q)-b)])
boost.dot(end-start)[-1].simplify()