Calculation of tail probabilities is of fundamental importance in several domains, such as for example risk assessment. One major challenge consists in the computation of low-failure probability and multiple-failure regions, especially when an unbiased estimation of the error is required. Methods developed in literature rely mostly on the construction of an adaptive surrogate, tackling some problems such as the metamodel building criterion and the global computational cost, at the price of a generally biased estimation of the failure probability. In this paper, we propose a novel algorithm permitting to both building an accurate metamodel and to provide a statistically consistent error. In fact, it relies on a novel metamodel building strategy, which aims to refine the limit-state region in all the branches "equally", even in the case of multiple failure regions, with a robust stopping building criterion.Secondly, two "quasi-optimal" importance sampling techniques are used, which permit, by exploiting the accurate knowledge of the metamodel, to provide an unbiased estimation of the failure probability, even if the metamodel is not fully accurate. As a consequence, the proposed method provides a very accurate unbiased estimation even for low failure probability or multiple failure regions.Several numerical examples are carried out, showing the very good performances of the proposed method with respect to the state-of-the-art in terms of accuracy and computational cost.Additionally, another importance sampling technique is proposed in this paper, permitting to drastically reduce the computational cost when estimating some reference values, or when a very weak failure-probability event should be computed directly from the metamodel.