You are using a very simplified treatment for the problem that ignores perhaps the most important part of it: the boundary conditions for the electromagnetic field in contact with the dielectric.
However, considering this treatment, your expression (5) for the Hamiltonian is incorrect. The Hamiltonian is not position dependent, and in your expression you expand it in momentum basis, so $\omega(p)$ cannot have different values depending on the position: $\omega$ does not depend on the position, otherwise the Hamiltonian would.
In this simplified formulation of the problem, the correct way of approaching this problem is by considering the Hamiltonian, given by
$$H = \int d z \left(\epsilon(z) :E(t,z)^2: + :B(t,z)^2:\right).$$
The normal ordering is necessary to prevent divergences, and is present even in the standard case, with no dielectric medium. Notice that $\epsilon(z)$ should be defned by \begin{equation} \epsilon(z) = \begin{cases}\epsilon, 0\leq z\leq L\ 1, \text{otherwise}\end{cases}. \end{equation}
After plugging in the values for $E$ and $B$ given by your expression (4), you will obtain a result that will look like
$$H = \frac{1}{4\pi}\int dp \:p\: \hat{a}^\dagger_p \hat{a}_p +(\epsilon-1)\int_0^L d z :E(t,z)^2:,$$
so that the modification to the dynamics of the system will come from the second term on the right hand side of the equation above. However, you would not be able to use expressions for the Dirac delta or other functions in that expression due to the finite interval, and it would end up looking quite bad: $$\begin{aligned} \int_0^L dz :E(t,z)^2: = \frac{1}{8\pi^2} \int_0^L dz \int dp \int dp'\Big(e^{-i(p+p')(t-z)}\hat{a}p\hat{a}{p'}+e^{-i(p-p')(t-z)}\hat{a}{p'}^\dagger\hat{a}{p}\+e^{i(p-p')(t-z)}\hat{a}p^\dagger\hat{a}{p'}+e^{i(p+p')(t-z)}\hat{a}p^\dagger\hat{a}{p'}^\dagger\Big). \end{aligned}$$
Using
$$\int_0^L dz e^{az} = \frac{1}{a}(e^{aL}-1)$$
in each one of the above integrals you would get your result. 526
(It looks a bit nicer making the cavity from $-L$ to $L$, so that everything becomes $\text{sinc}$ functions). In essence, there is no way of handling the quadratic terms on the creation and annihilation operators, which will make it impossible to obtain a simple solution for your state. Perturbation theory could be an option for finite times...
In any case, the appropriate way of handling this problem is to impose boundary conditions for the electric field inside the cavity. Under these conditions the field modes inside it would be discrete, while the modes outside would be continuous. However, it is very important to notice that the notion of photon wavepacket you have defined there depends explicitly on whether your field is within a cavity, or in infinite space. This is a general feature of particle excitations in quantum field theory: they are not local, and depend on the mode decomposition chosen for the field, which depends on the boundary conditions. In fact, defining particle states for the case of fields with general boundary conditions is not an easy problem to tackle. 795