For a toric pair $(X, D)$, where $X$ is a projective toric variety of dimension $d-1geq 1$ and $D$ is a very ample $T$-Cartier divisor, we show that the Hilbert-Kunz density function $HKd(X, D)(lambda)$ is the $d-1$ dimensional volume of ${overline {mathcal P}}_D cap {z= lambda}$, where ${overline {mathcal P}}_Dsubset {mathbb R}^d$ is a compact $d$-dimensional set (which is a finite union of convex polytopes). We also show that, for $kgeq 1$, the function $HKd(X, kD)$ can be replaced by another compactly supported continuous function $varphi_{kD}$ which is `linear in $k$. This gives the formula for the associated coordinate ring $(R, {bf m})$: $$lim_{kto infty}frac{e_{HK}(R, {bf m}^k) - e_0(R, {bf m}^k)/d!}{k^{d-1}} = frac{e_0(R, {bf m})}{(d-1)!}int_0^inftyvarphi_D(lambda)dlambda, $$ where $varphi_D$ (see Proposition~1.2) is solely determined by the shape of the polytope $P_D$, associated to the toric pair $(X, D)$. Moreover $varphi_D$ is a multiplicative function for Segre products. This yields explicit computation of $varphi_D$ (and hence the limit), for smooth Fano toric surfaces with respect to anticanonical divisor. In general, due to this formulation in terms of the polytope $P_D$, one can explicitly compute the limit for two dimensional toric pairs and their Segre products. We further show that (Theorem~6.3) the renormailzed limit takes the minimum value if and only if the polytope $P_D$ tiles the space $M_{mathbb R} = {mathbb R}^{d-1}$ (with the lattice $M = {mathbb Z}^{d-1}$). As a consequence, one gets an algebraic formulation of the tiling property of any rational convex polytope.