We perform a multivariate analysis of Higgs-pair production in HH→bb¯γγ channel at the HL-LHC to probe the trilinear Higgs self-coupling (THSC) λ3H, which takes the value of 1 in the standard model (SM). We consider all known background processes. For the signal we are the first to adopt the most recent event generator of POWHEG-BOX-V2 to exploit the next-to-leading order (NLO) distributions for the toolkit for multivariate data analysis (TMVA), taking account of the full top-quark mass dependence. Through boosted decision tree (BDT) analysis trained for λ3H=1, we find that the significance can reach up to 1.95 with about 9 signal and 18 background events. In addition, the THSC can be constrained to 1.00<λ3H<6.22 at 95% confidence level (CL). We also perform a likelihood fitting of Mγγbb distribution and find the 1σ confidence interval (CI) of 0.1<λ3H<2.2∪5.4<λ3H<6.6 for the λ3H=1 nominal set. On the other hand, using BDTs trained for each value of λ3H, we find a bulk region of 0.5≲λ3H≲4.5, for which it is hard to pin down the THSC.