  We perform a multivariate analysis of Higgs-pair production in $HH \to b\bar
b \gamma\gamma$ channel at the HL-LHC to probe the trilinear Higgs
self--coupling $\lambda_{3H}$, which takes the value of 1 in the SM. We
consider all the known background processes. Also, for the signal we are the
first to adopt the most recent event generator of POWHEG-BOX-V2 to exploit the
NLO distributions for Toolkit for Multivariate Data Analysis (TMVA), taking
account of the full top--quark mass dependence. Through Boosted Decision Tree
(BDT) analysis trained for $\lambda_{3H}=1$, we find that the significance can
reach up to 1.95 with about $9$ signal and $18$ background events. In addition,
the Higgs boson self-coupling can be constrained to $1.00 < \lambda_{3H} <
6.22$ at 95\% confidence level (CL). We also perform a likelihood fitting of
$M_{\gamma\gamma bb}$ distribution and find the $1\sigma$ confidence interval
(CI) of $0.1 < \lambda_{3H} < 2.2\, \cup\,5.4 < \lambda_{3H} < 6.6$ for the
$\lambda_{3H}=1$ nominal set. On the other hand, using BDTs trained for each
value of $\lambda_{3H}$, we find a bulk region of $0.5 <\sim \lambda_{3H}
<\sim~ 4.5$, for which it is hard to pin down the trilinear coupling.
