The precise estimation of the mass of galaxy clusters is a major issue for cosmology. Large galaxy cluster surveys rely on scaling laws that relate cluster observables to their masses. From the high-resolution observations of ∼45 galaxy clusters with the NIKA2 and XMM-Newton instruments, the NIKA2 Sunyaev-Zel’dovich Large Program should provide an accurate scaling relation between the thermal Sunyaev-Zel’dovich effect and the hydrostatic mass. In this paper we present an exhaustive analysis of the hydrostatic mass of the well-known galaxy cluster CL J1226.9+3332, the highest-redshift cluster in the NIKA2 Sunyaev-Zel’dovich Large Program at z = 0.89. We combined the NIKA2 observations with thermal Sunyaev-Zel’dovich data from the NIKA, Bolocam, and MUSTANG instruments and XMM-Newton X-ray observations, and tested the impact of the systematic effects on the mass reconstruction. We conclude that slight differences in the shape of the mass profile can be crucial when defining the integrated mass at R₅₀₀, which demonstrates the importance of the modelling in the mass determination. We prove the robustness of our hydrostatic mass estimates by showing the agreement with all the results found in the literature. Another key factor for cosmology is the bias of the masses estimated assuming the hydrostatic equilibrium hypothesis. Based on the lensing convergence maps from the Cluster Lensing And Supernova survey with Hubble (CLASH) data, we obtain the lensing mass estimate for CL J1226.9+3332. From this we are able to measure the hydrostatic-to-lensing mass bias for this cluster, which spans from 1 − b_HSE/lens ∼ 0.7 to 1, presenting the impact of data sets and mass reconstruction models on the bias.