For $c=0$ this result is knows as reflection principle (see e.g. René Schilling/Lothar Partzsch: Brownian Motion - An Introduction to Stochastic Processes, Chapter 6) and follows from the Markov property and symmetry of Brownian motion. However, for $c>0$ the proof is more involved since we have to get rid of the drift term.
Since by definition
$$[H_a \leq t] = \left[ \sup_{s \leq t} X_s \geq a \right] \tag{1}$$
determining the distribution of $H_a$ is equivalent to finding the distribution of $\sup_{s \leq t} X_s$. In order to find the distribution of the latter, we need two ingredients: Girsanov's theorem and the joint distribution $(X_t,\sup_{s \leq t} X_s)$ for a Brownian motion $(X_t)_{t \ge 0}$.
Girsanov theorem: Let $c \in \mathbb{R}$ and $(B_t)_{t \geq 0}$ be a Brownian motion on a probability space $(\Omega,\mathcal{A},\mathbb{P})$. Then $$X_t := B_t+ct, \qquad t \leq T,$$ is a Brownian motion with respect to the probability measure $$d\mathbb{Q} := d\mathbb{Q}_T := \exp \left( -c B_T - \frac{c^2}{2} T \right) d\mathbb{P}.$$
For a proof see e.g. René Schilling/Lothar Partzsch: Brownian Motion - An Introduction to Stochastic Processes, Chapter 18.
Joint distribution of $(X_t, \sup_{s \leq t} X_s)$: Let $(X_t)_{t \geq 0}$ be a Brownian motion on a probability space $(\Omega,\mathcal{A},\mathbb{Q})$. Then the joint distribution $(X_t,\sup_{s \leq t} X_s)$ equals $$\mathbb{Q} \left[ X_t \in dx, \sup_{s \leq t} X_s \in dy \right] = \frac{2 (2y-x)}{\sqrt{2\pi t^3}} \exp \left(- \frac{(2y-x)^2}{2t} \right) 1_{[-\infty,y]}(x) \, dx \, dy. \tag{2}$$
For a proof see e.g. René Schilling/Lothar Partzsch: Brownian Motion - An Introduction to Stochastic Processes, Exercise 6.8 (there are full solutions available on the web).
So let's finally put it all together: It follows from $(1)$ and the definition of the probability measure $\mathbb{Q}_T$ that
$$\begin{align*} \mathbb{P}[H_a \leq T] &= \mathbb{P} \left[ \sup_{s \leq T} X_s \geq a \right] = \int 1_{[a,\infty)} \left( \sup_{s \leq T} X_s \right) \, d\mathbb{P} \\ &= \int 1_{[a,\infty)}\left( \sup_{s \leq T} X_s \right) \exp \left( c B_T + \frac{c^2}{2} T \right) \, d\mathbb{Q}_T \\ &= \int 1_{[a,\infty)}\left( \sup_{s \leq T} X_s \right)\exp\left(c X_T- \frac{c^2}{2} T \right) \, d\mathbb{Q}_T. \end{align*}$$
By Girsanov's theorem, $(X_t)_{t \leq T}$ is a Brownian motion with respect to $\mathbb{Q}_T$ and therefore $(2)$ gives
$$\begin{align*} \mathbb{P}[H_a \leq T] &= \exp\left(- \frac{c^2}{2} T \right) \int_{y \geq a} \int_{x \leq y} e^{cx} \frac{2 (2y-x)}{\sqrt{2\pi T^3}} \exp \left(- \frac{(2y-x)^2}{2T} \right) \, dx \, dy. \end{align*}$$
It remains to calculate the integral expression. First of all, by Fubini's theorem,
$$\begin{align*} \mathbb{P}[H_a \leq T] &= \frac{1}{\sqrt{2\pi T}} \exp \left(- \frac{c^2}{2} T \right) \left(\int_{x \geq a} e^{cx} I_1(x) \, dx + \int_{x \leq a} e^{-cx} I_2(x) \, dx \right) \\ &:=J_1+J_2 \tag{3} \end{align*}$$
where
$$\begin{align*} I_1(x):= \int_{y \geq x} \frac{2(2y-x)}{T} \exp \left(- \frac{(2y-x)^2}{2T} \right) \, dy &= \left[ - \exp \left(- \frac{(2y-x)^2}{2T} \right) \right]_{y=x}^{\infty} \\ &= \exp \left(- \frac{x^2}{2T} \right) \\ I_2(x) := \int_{y \geq a} \frac{2(2y-x)}{T} \exp \left(- \frac{(2y-x)^2}{2T} \right) \, dy &=\exp \left(- \frac{(2a-x)^2}{2T} \right). \end{align*}$$
Hence,
$$\begin{align*} J_1 &= \frac{1}{\sqrt{2\pi T}} \exp \left(- \frac{c^2}{2} T \right) \int_{x \geq a} e^{cx} I_1(x) \, dx \\ &= \frac{1}{\sqrt{2\pi T}} \int_{x \geq a} \exp \left(- \frac{(x-cT)^2}{2T} \right) \, dx \\ &= \frac{1}{\sqrt{\pi}} \int_{z \geq \frac{a-cT}{\sqrt{2T}}} \exp(-z^2) \, dz \tag{4} \end{align*}$$
and
$$\begin{align*} J_2 &= \frac{1}{\sqrt{2\pi T}} \exp \left(- \frac{c^2}{2} T \right) \int_{x \leq a} e^{cx} I_2(x) \, dx \\ &\stackrel{u:=2a-x}{=} \frac{1}{\sqrt{2\pi T}} \exp \left(- \frac{c^2}{2} T \right) \int_{u \geq a} e^{c(2a-u)} \exp\left(-\frac{u^2}{2T} \right) \, du \\ &=\ldots = \frac{e^{2ac}}{\sqrt{\pi}} \int_{z \geq \frac{a+CT}{\sqrt{2T}}} \exp(-z^2) \, dz. \tag{5} \end{align*}$$
Now if we differentiate $(3)$ with respect to $T$, using $(4)$ and $(5)$, we get
$$\begin{align*} \frac{d}{dT} \mathbb{P}(H_a \leq T) &= \frac{-1}{\sqrt{\pi}} \exp \left( - \frac{(a-cT)^2}{2T} \right) \left( \frac{-c}{\sqrt{2T}} - \frac{a-cT}{2 \sqrt{2} T^{3/2}} \right) \\ &\quad + \frac{-1}{\sqrt{\pi}} \underbrace{e^{2ac} \exp \left( - \frac{(a+cT)^2}{2T} \right)}_{\exp(-(a-cT)^2/2T)} \left( \frac{c}{\sqrt{2T}} - \frac{a+cT}{2 \sqrt{2} T^{3/2}} \right) \\ &= \frac{a}{\sqrt{2\pi T^3}} \exp \left(- \frac{(a-cT)^2}{2T} \right). \end{align*}$$
There's actually two definitions of solution to SDE. Strong and weak.
Strong solution. Given a probability space $
(\Omega, \mathcal F, \mathcal F_t,P)$ and a Brownian motion $B(t)$ on that space adapted to $\mathcal F_t$, a solution to $dX=fdt+g~dB(t)$ with initial condition $x\in \mathbb R$ is a continuous process adapted to $\mathcal F_t$ such that $X(t)=x+\int_0^t fdt+\int_0^t g~dB(t)$.
There is also
Weak solution. A process $X(t)$ is called a weak solution to $dX=fdt+g~dB(t)$ if there exists a probability space $
(\Omega, \mathcal F, \mathcal F_t,P)$ and a Brownian motion $B(t)$ on that space adapted to $\mathcal F_t$ such that $X(t)=x+\int_0^t fdt+\int_0^t g~dB(t)$.
Note that strong solution implies weak solution.
As an example of something that has weak but not strong solution, consider Tanaka's equation
$$dX(t)=\text{sign}(X(t))dB(t)$$
It can be shown that this has no strong solution by considering local times. However let $X(t)$ be a Brownian motion and then $X(t)$ is a weak solution. Note that $\int_0^t\text{sign}^2(s)ds=t<\infty$ so the Ito integral $\int_0^t\frac{1}{\text{sign}(s)}dX(s)=\int_0^t \text{sign}(s) dX(s)$ exists. Noting that the quadratic variation of this integral is $t$, and Ito integrals are continuous martingales, so we know $\int_0^t\frac{1}{\text{sign}(s)}dX(s)=\tilde{B}(t)$ is a Brownian motion.
Thus, $dX(t)=\text{sign}(t)d\tilde{B}(t)$.
The difference here is that the Brownian motion is constructed in terms of our solution not the other way around.
Best Answer
If $(B_t^x)_{t \geq 0}$ is a Brownian motion started at $x \in \mathbb{R}^d$, then the density of $B_t^x$ is given by
$$p_t^x(y) = \frac{1}{(2\pi t)^{d/2}} \exp \left(- \frac{|y-x|^2}{2t} \right).$$
Using the Dirac measure we can rewrite this as follows:
$$p_t^x(y) = \frac{1}{(2\pi t)^{d/2}} \int \exp \left(- \frac{|y-z|^2}{2t} \right) \delta_x(dz) $$
Note that $\delta_x$ is the initial distribution of the Brownian motion $(B_t^x)_{t \geq 0}$. If we replace $\delta_x$ by some general distribution, say $\mu$, we get
$$p_t^{\mu}(y) = \frac{1}{(2\pi t)^{d/2}} \int \exp \left(- \frac{|y-z|^2}{2t} \right) \, \mu(dz).$$
Roughly this means that we mix the densities $(p_t^x)_{x \in \mathbb{R}^d}$ according to our given initial distribution. In particular, if $\mu(dz) = u_0(z) \, dz$ for some density $u_0$, we have
$$p_t^{\mu}(y) = \frac{1}{(2\pi t)^{d/2}} \int \exp \left(- \frac{|y-z|^2}{2t} \right) u_0(z) \, dz.$$
A straight-forward calculation shows that this function does indeed satisfy the heat equation. Moreover, one can show that $p_t^{\mu}(y) \to u_0(y)$ as $t \to 0$. Note that the density $(p_t^{\mu})$ is the convolution of the initial distribution and the heat kernel.
This is the probabilistic approach. The more popular approach to solve the heat equation uses Fourier methods (i.e. take the Fourier transform of the heat equation, do some calculations and then invert the Fourier transform to get the solution $u$). You will find this in (almost) any book on PDEs.