  In a cylinder $\Omega_T=\Omega\times (0,T)\subset \R^{n+1}_+$ we study the
boundary behavior of nonnegative solutions of second order parabolic equations
of the form \[ Hu =\sum_{i,j=1}^ma_{ij}(x,t) X_iX_ju - \p_tu = 0, \
(x,t)\in\R^{n+1}_+, \] where $X=\{X_1,...,X_m\}$ is a system of $C^\infty$
vector fields in $\Rn$ satisfying H\"ormander's finite rank condition
\eqref{frc}, and $\Omega$ is a non-tangentially accessible domain with respect
to the Carnot-Carath\'eodory distance $d$ induced by $X$. Concerning the
matrix-valued function $A=\{a_{ij}\}$, we assume that it be real, symmetric and
uniformly positive definite. Furthermore, we suppose that its entries $a_{ij}$
be H\"older continuous with respect to the parabolic distance associated with
$d$. Our main results are: 1) a backward Harnack inequality for nonnegative
solutions vanishing on the lateral boundary (Theorem \ref{T:back}); 2) the
H\"older continuity up to the boundary of the quotient of two nonnegative
solutions which vanish continuously on a portion of the lateral boundary
(Theorem \ref{T:quotients}); 3) the doubling property for the parabolic measure
associated with the operator $H$ (Theorem \ref{T:doubling}). These results
generalize to the subelliptic setting of the present paper, those in Lipschitz
cylinders by Fabes, Safonov and Yuan in [FSY] and [SY]. With one proviso: in
those papers the authors assume that the coefficients $a_{ij}$ be only bounded
and measurable, whereas we assume H\"older continuity with respect to the
intrinsic parabolic distance.
