Neutral beam injection (NBI) will be one of the main sources of heating and non-inductive current drive in ITER. Due to high level of injected power the beam induced heat loads can, however, jeopardize the integrity of the first wall of the device, particularly in the presence of non-axisymmetric perturbations of the magnetic field. Neutral beam injection can also destabilize Alfvn eigenmodes and energetic particle modes (EPM), and act as a source of plasma rotation. Therefore, reliable and accurate simulation of NBI is important formaking predictions for ITER as well as any other current or future fusion device. This paper introduces a new beamletbased neutral beam ionization model called BBNBI. It takes into account the fine structure of the injector, follows the injected neutrals until ionization, and generates a source ensemble of ionized NBI test particles for slowing down calculations. BBNBI can be used as a stand-alone model but together with the particle following code ASCOT, it forms a complete and sophisticated tool for simulating neutral beam injection. The test particle ensembles fromBBNBI are found to agree well with those produced by PENCIL for JET, and those produced by NUBEAM both for JET and ASDEX Upgrade plasmas. The first comprehensive comparisons of beam slowing down profiles of interest from BBNBI+ASCOT with results from PENCIL and NUBEAM/TRANSP, for both JET and AUG, are presented. It is shown that, for an axisymmetric plasma, BBNBI+ASCOT and NUBEAM agree remarkably well. Together with earlier 3D studies, these results further validate using BBNBI+ASCOT also for studying phenomena that require particle following in a truly three-dimensional geometry.