  Fragmentation functions (FFs) are essential non-perturbative QCD inputs for
predicting hadron production cross sections in high energy scatterings. In this
study, we present a joint determination of FFs for light charged hadrons
through a global analysis at next-to-leading order (NLO) in QCD. Our analysis
incorporates a wide range of precision measurements from the LHC, as well as
data from electron-positron collisions and semi-inclusive deep inelastic
scatterings. By including measurements of jet fragmentation at the LHC in our
global analysis, we are able to impose strong constraints on the gluon FFs. A
careful selection of hadron kinematics is applied to ensure the validity of
factorization and perturbative calculations of QCD. In addition, we introduce
several methodological advances in fitting, resulting in a flexible
parametrization form and the inclusion of theoretical uncertainties from
perturbative calculations. Our best-fit predictions show very good agreement
with the global data, with $\chi^2/N_{pt}\sim 0.90$. We also generate a large
number of Hessian error sets to estimate uncertainties and correlations of the
extracted FFs. FFs to charged pions (kaons and protons) are well constrained
for momentum fractions down to 0.01 (0.1). Total momentum of partons carried by
light charged hadrons are determined precisely. Their values for $u$, $d$
quarks and gluon saturate at about 50\% for a lower cut of the momentum
fraction of 0.01. Pulls from individual datasets and impact of various choices
of the analysis are also studied in details. Additionally, we present an update
of the FMNLO program used for calculating hadron production cross sections. Our
FFs, including the error sets (denoted as NPC23), are publicly available in the
form of LHAPDF6 grids.
