First of all, I would like to thank John Dawkins for his kind consideration and time. Honestly, I could not convince myself by reading the sources John preferred me. So, I decided to search more and find a solution to my question. Here, is my finding.
At the beginning, let's consider a standard Brownian motion $B(t)$. Assume that $T_a = inf\{t: B(t) = a\}$ stands for hitting time, i.e., the first time when Brownian motion hits level $a>0$, and $M(t) = \sup_{0\leqslant s\leqslant t}B(s)$. We have
\begin{eqnarray}
P(B(t)> a) = P(B(t)> a , M(t)> a)+P(B(t)> a , M(t)\leqslant a)
\end{eqnarray}
Note that $P(B(t)> a , M(t)\leqslant a)$ since $M(t)\geqslant B(t)$. Now
\begin{equation*}
P(B(t)> a , M(t)> a)=
P(B(t)> a | M(t)> a)\times P(M(t)> a)\\ =
P(B(t)> a | T_a< t)\times P(M(t)> a)
\end{equation*}
we have that $B(T_a+s)-a$ is a Brownian motion. Conditioned on $ T_a\leqslant t$,
\begin{eqnarray}
P(B(t)> a | T_a< t) = P(B(T_a+(t-T_a))-a> 0 |T_a< t) = \frac{1}{2}
\end{eqnarray}
since the Brownian motion satisfies $P(B(t)> 0) = \frac{1}{2}$ for every t. Applying this identity, we obtain
\begin{eqnarray}
P(B(t)> a) = \frac{1}{2}P(M(t)> a)
\end{eqnarray}
or
\begin{eqnarray}
P(M(t)> a) = 2P(B(t)>a) = P(|B(t)|>a) \Rightarrow P(M(t)\leqslant a)= P(|B(t)|\leqslant a)
\end{eqnarray}
we know that $B(t)\sim N(0,t)$, so
\begin{eqnarray*}
P(M(t)\leqslant a) =
\frac{1}{\sqrt{2\pi t}}\int_{-a}^{a} e^{-\frac{x^2}{2t}}
= \frac{2}{\sqrt{2\pi t}}\int_{0}^{a} e^{-\frac{x^2}{2t}}
\end{eqnarray*}
we can also write it down in this way
\begin{eqnarray}
P(M(t)\leqslant a) = 2P(B(t)\leqslant a) = 2\Phi(\frac{a}{\sqrt{t}})
\end{eqnarray}
Now, consider a stock price $S_t = S_0 e^{\sigma W_t}$,i.e., a Brownian motion process without drift $(\mu = 0)$. Then, we conclude that
\begin{eqnarray}
M_0(T) = \max_{0\leqslant t \leqslant T}S_t = \max_{t\in[0,T]}S_0 e^{\sigma W_t} = S_0e^{\sigma\max_{t\in[0,T]}W_t} = S_0e^{\sigma M(T)}
\end{eqnarray}
Therefore,
\begin{eqnarray*}
P(M_0(T)\leqslant y)&=&P(S_0e^{\sigma M(T)}\leqslant y) \\
&=& P(e^{\sigma M(T)}\leqslant \frac{y}{S_0}) \\
&=& P(\sigma M(T)\leqslant log(\frac{y}{S_0})) \\
&=& 2\Phi\left(\frac{log(\frac{y}{S_0})}{\sqrt{T}}\right),\quad y>S_0
\end{eqnarray*}
and for the inverse distribution function, we get this result
\begin{eqnarray}
F^{-1}_{M_0(T)}(x) = S_0e^{\sqrt{T}\Phi^{-1}(\frac{x}{2})}
\end{eqnarray}
Generally, for the geometric Brownian motion, it is more complicated to find the exact formula. I only put the final result. Considering the following Geometric Brownian motion for a stock price,
\begin{eqnarray*}
S(t) = S_0exp\left((r-d-\frac{\sigma^2}{2})t+\sigma W(t)\right)
\end{eqnarray*}
and maximum process of such model
\begin{eqnarray}
M_0(T) = \max_{0 \leqslant t \leqslant T}S(t)
\end{eqnarray}
we have that
\begin{eqnarray}
P(M_0(T)\leqslant x) = \Phi\left(\frac{-(r-d-\frac{\sigma^2}{2})T+log(\frac{x}{S_0})}{\sigma \sqrt{T}}\right)-\left(\frac{S_0}{x}\right)^{1-2\frac{(r-d)}{\sigma^2}}\Phi\left(\frac{-(r-d-\frac{\sigma^2}{2})T-log(\frac{x}{S_0})}{\sigma \sqrt{T}}\right), \quad x>S_0
\end{eqnarray}
