Using hadron-in-jet data in a global analysis of D^{*} fragmentation functions

# Using hadron-in-jet data in a global analysis of D∗ fragmentation functions

## Abstract

We present a novel global QCD analysis of charged -meson fragmentation functions at next-to-leading order accuracy. This is achieved by making use of the available data for single-inclusive -meson production in electron-positron annihilation, hadron-hadron collisions, and, for the first time, in-jet fragmentation in proton-proton scattering. It is shown how to include all relevant processes efficiently and without approximations within the Mellin moment technique, specifically for the in-jet fragmentation cross section. The presented technical framework is generic and can be straightforwardly applied to future analyses of fragmentation functions for other hadron species, as soon as more in-jet fragmentation data become available. We choose to work within the Zero Mass Variable Flavor Number Scheme which is applicable for sufficiently high energies and transverse momenta. The obtained optimum set of parton-to- fragmentation functions is accompanied by Hessian uncertainty sets which allow one to propagate hadronization uncertainties to other processes of interest.

## I Introduction

Cross sections at collider experiments can often be reliably calculated within the framework of perturbative Quantum Chromodynamics (pQCD). The crucial foundation for such computations are so-called factorization theorems that allow for a systematic separation of perturbative and non-perturbative physics ref:fact (). Well-known examples for the latter are parton distribution functions (PDFs) that are, by now, rather tightly constrained by global QCD fits to data and are a crucial asset in all scattering processes with hadrons in the initial-state.

Whenever an observable involves detected hadrons in the final-state, the theoretical calculation requires another type of non-perturbative functions as input. These parton-to-hadron fragmentation functions (FFs) describe the non-perturbative transition of a parton produced in the hard-scattering event into the observed hadron. Like PDFs, these functions were shown to be universal and can be only extracted from data through global QCD analyses. The knowledge of FFs for different hadron species and estimates of their uncertainties is therefore vital for precise theoretical calculations and, hence, has received quite some interest in the past; see, for instance, Ref. Metz:2016swz () for a recent review.

In this work, we consider the hadronization of quarks and gluons into heavy-flavored mesons, more specifically, charged -mesons, that are of particular relevance in the era of the LHC. In general, the theoretical treatment of heavy quarks itself provides a unique laboratory to test pQCD. Correctly describing heavy flavor cross sections that have been measured both at very high energies at the LHC and at various low energy experiments poses unique challenges to our understanding of QCD. Charm production cross sections are used, for example, to constrain the gluon PDF at small- Gauld:2015yia (), and they play a vital role in cosmic-ray and neutrino astrophysics ref:astro (). Another important area of research concerns the modification of heavy flavor yields in heavy-ion collisions ref:heavyion () where highly energetic partons can traverse the quark-gluon plasma thereby attaining valuable information about the properties of the QCD medium. For instance, the energy loss mechanisms, that allow for a quantitative description of in-medium effects, crucially depend on the underlying fragmentation process.

In pQCD calculations, the heavy quark mass introduces an additional large scale apart from some other hard scale characterizing the process, such as a measured transverse momentum . These multi-scale problems carry additional theoretical challenges as compared to processes involving only light quarks and gluons. There are various approaches in the literature of how to deal with heavy quark masses in general and in the fragmentation process in particular. In the context of collisions relevant for LHC phenomenology the following schemes have been put forward and used in their various kinematic regimes of applicability. In the Fixed Flavor Number Scheme (FFNS) ref:ffns (), the heavy quark is not treated as an active parton in the proton but, instead, is solely produced extrinsically in the hard scattering. Logarithms of the ratio of the heavy quark mass and the hard scale of the process, , are only taken into account in fixed order perturbation theory. Therefore, this scheme is applicable in the region . The Zero Mass Variable Flavor Number Scheme (ZMVFNS), on the other hand, is only applicable in the limit . Here, large logarithms of are resummed through DGLAP evolution equations to all orders. is set to zero in all partonic cross sections, and the heavy quark is treated as an active, massless parton in the proton.

In the context of fragmentation processes, the Fixed Order plus Next-to-Leading Logarithms prescription (FONLL) ref:cacciari (); Cacciari:2005uk () as well as the General Mass Variable Flavor Number Scheme (GMVFNS) ref:gmvfns (); Kneesch:2007ey () are examples of unified frameworks to cover both the high region, , and the low tail, , similar to the ZMVFNS and FFNS, respectively. In the FONLL approach, the FFs of heavy-flavored mesons are separated into a perturbatively calculable parton-to-heavy quark contribution and a non-perturbative heavy quark-to-heavy meson piece that is fitted to data. This separation is possible as the heavy quark mass sets an additional scale in the perturbative regime. Instead, in the GMVFNS the entire parton-to-heavy meson FF is treated as a non-perturbative function and is extracted from the available data. We note that another scheme was developed recently in Ref. Fickinger:2016rfd () within the framework of Soft Collinear Effective Theory (SCET).

Since we are primarily interested in LHC phenomenology in this work, in particular the impact of in-jet fragmentation data at , we choose to work in the ZMVFNS using purely non-perturbative FFs similar to the analyses of FFs for light hadron species. As will be discussed in detail below, the inclusive -spectrum of charged -mesons in collisions can be fairly well described in the ZMVFNS down to rather low values of about in spite of imposing a cut when fitting data.

Traditionally, the main reference process to determine FFs is semi-inclusive electron-positron annihilation (SIA), . Here, denotes the detected hadron and the unobserved final-state remnant. To the best of our knowledge, all the approaches to heavy quark fragmentation mentioned above rely only on SIA data to determine the relevant non-perturbative input following similar non-global fits of light hadron (pion, kaon) FFs ref:lightffs (). While quark-to-hadron FFs can be relatively well constrained from SIA data, it is, in particular, the gluon-to-hadron FF that is at best only very poorly constrained by SIA data alone. Therefore, global QCD analyses of light hadron FFs have also included vital proton-proton scattering data, , in order to better constrain the gluon FF. In addition, Semi-Inclusive Deep Inelastic Scattering (SIDIS), , data are needed to perform a quark-antiquark and quark flavor separation of FFs. Such global fits of light hadron FFs can be found in ref:dss ().

In this paper, we will provide the first global QCD analysis of charmed-meson FFs following the framework outlined by the DSS group in ref:dss () at next-to-leading order (NLO) accuracy using the Mellin moment technique Stratmann:2001pb (). We note that recently first efforts have been made to perform fits of light hadron FFs at next-to-next-to-leading order (NNLO) accuracy Anderle:2015lqa (); Bertone:2017tyb (), and also by including all-order resummations Anderle:2012rq (); Anderle:2016czy (). So far, these efforts have been limited to SIA data only due to the lack of other single-inclusive particle production cross sections at NNLO accuracy; see, for example, Anderle:2016kwa () for the progress of an ongoing SIDIS calculation at NNLO. As has become customary for both PDF and FF analyses these days, we also present an attempt to estimate the remaining uncertainties of the extracted FF, for which we adopt the Hessian method ref:hessian (). The Hessian uncertainty sets can be used to propagate hadronization uncertainties to any other processes of interest such as, for instance, high- -meson production in proton-nucleus collisions at the LHC or BNL-RHIC.

Besides the processes that are traditionally included in global analyses of FFs, like SIA and inclusive high- hadron production in collisions, we also include for the first time in-jet fragmentation data from the LHC. Specifically, we include data for the “jet fragmentation function”, , where a hadron is identified inside a jet. We consider the observable, where the longitudinal momentum distribution differential in is measured, with () denoting the hadron (jet) transverse momentum. The fact that at leading order (LO) accuracy the in-jet observable is directly probing the dependence of FFs explains their potential relevance for analyses of FFs. In-jet fragmentation was pioneered in ref:injet () for exclusive jet samples. The extension to inclusive jet samples was developed in Arleo:2013tya (); Kaufmann:2015hma () within standard pQCD at NLO accuracy, allowing for a direct comparison with data from the LHC. In Ref. Kang:2016ehg () the result was re-derived within SCET. Thanks to the effective field theory treatment, the additional all-order resummation of single logarithms in the jet-size parameter was achieved, yielding consistent results at NLO+NLL accuracy. In this work, we will work at NLO accuracy, and we leave a detailed study of the impact of NLO+NLL corrections on fits of FFs for a dedicated, future publication.

In Ref. Chien:2015ctp (), it was found that the -in-jet data from ATLAS Aad:2011td () are not well described by existing fits of -meson FFs Kneesch:2007ey () even though they give a good description of both SIA and inclusive data; see Ref. Bain:2016clc () for related work. This leads to the important question, which we address in detail in this work, if there is a real tension between the fitted data sets and the in-jet observable or if it is possible to accommodate all data sets in a combined, global fit. We note that apart from the -in-jet data by ATLAS Aad:2011td () there are also in-jet results available from the LHC for unidentified light charged hadrons ref:unid-data (), mainly in heavy-ion collisions though, as well as for prompt and non-prompt production in jets Aaij:2017fak (). In this first exploratory study of the impact of in-jet data on fits of FFs we therefore limit ourselves to developing the necessary theoretical framework and to a global analysis of parton-to- and parton-to- FFs utilizing the ATLAS in-jet data. However, we wish to emphasize that the technical framework presented below is generic and can be straightforwardly applied to future analyses of fragmentation functions for other hadron species as soon as more in-jet fragmentation data become available.

Finally, we notice that various combined differential cross section data for charged mesons obtained in deep-inelastic lepton-proton collisions are available from the H1 and ZEUS Collaborations H1:2015dma (). Since the data extend down to relatively low values of transverse momentum and photon virtuality , they need to be described in a theoretical framework which keeps the full dependence on the charm quark mass H1:2015dma (). Hence, these data cannot be included in our current global QCD analysis that is based on the ZMVFNS approximation.

The remainder of this paper is organized as follows. In Section II, we discuss the technical framework for all three processes that are included in our global QCD analysis, namely , , and , with particular emphasis on the latter. Note that throughout this paper, collectively denotes both charged mesons, i.e.  and/or . Next, in Section III, we briefly present the details of our analysis comprising the parametrization of the FFs at some input scale, the selection of experimental data and cuts imposed on the fit, the Mellin moment technique used throughout this paper, and the Hessian uncertainty method. In Section IV, we present and discuss the results of our global analysis of parton-to- FFs at NLO accuracy, and compare the results of the fit to the available data. In addition, we compare to the previous fit provided by Ref. Kneesch:2007ey (). In Section V, we draw our conclusions and present a brief outlook.

## Ii Technical Framework

### ii.1 Single-inclusive e+e− annihilation

The cross section for the single-inclusive annihilation process, , is usually normalized to the total hadronic cross section and may be written schematically as

 1σtotdσe+e−→hXdz=σ0σtot[FhT(z,Q2)+FhL(z,Q2)]. (1)

It is common to decompose the cross section (1) into a transverse () and longitudinal () part although this is of no practical relevance for -meson production. We have introduced the scaling variable

 z≡2Ph⋅qQ2\lx@stackrelc.m.s.=2EhQ, (2)

where and are the four momenta of the observed hadron and time-like boson, respectively. Moreover, . As is indicated in Eq. (2), reduces to the hadron’s energy fraction in the center-of-mass system (c.m.s.) frame and is often also labeled as ref:nason (). The total cross section for hadrons at NLO accuracy reads

 σtot=∑q^e2qσ0[1+αs(Q2)π], (3)

where and are the electromagnetic and the strong coupling, respectively, and . We denote the electroweak quark charges by , which may be found, for instance, in App. A of Ref. deFlorian:1997zj ().

To make factorization explicit, the transverse and longitudinal time-like structure functions in Eq. (1) can be written as a convolution of perturbative coefficient functions ,  ref:sia-coefficients (), and non-perturbative FFs ,

 Fhk(z,Q2) = ∑q^e2q{[Ckq⊗(Dhq+Dh¯q)](z,Q2) (4)

where . The standard convolution integral with respect to the first argument is denoted by the symbol and reads

 [f⊗g](z,…)≡∫10dx∫10dyf(x,…)g(y,…)δ(z−xy) (5)

for two arbitrary functions and .

As always, the notion of factorization as applied in Eq. (4) is only valid up to corrections proportional to inverse powers of the hard scale ref:nason (). For a one-scale process like SIA, the hard scale should be chosen to be of and itself should be at least of (few GeV). For simplicity, we have chosen the factorization and renormalization scales in Eq. (4) equal to the hard scale, i.e., .

Kinematical effects related to the non-zero mass of the produced hadron are another source of corrections to the factorized framework where is neglected throughout. Deviations of the data from theory are expected to show up at the lower end of the -spectrum, as we shall see in the phenomenological section, and are more pronounced for heavier than for light mesons. One usually introduces a cut in global analyses of FFs ref:lightffs (); ref:dss () below which the data cannot be used and the theory outlined above is not valid. Such a cut also avoids the region in where fixed-order evolution kernels receive large logarithmic corrections which otherwise can only be dealt with by all-order resummations, see, for instance, Ref. Anderle:2016czy ().

### ii.2 Single-inclusive D∗ production pp collisions

The production of high- hadrons in hadronic collisions offers valuable and complementary information compared to SIA data in global QCD analyses of FFs. The dominance of the and partonic subprocesses at not too large values of gives access to the gluon-to-hadron fragmentation function, which is only very indirectly accessible in SIA through scaling violations and, hence, largely unconstrained.

In addition to data for the sum of charged mesons, , from ATLAS Aad:2015zix () and LHCb Aaij:2013mga (); Aaij:2015bpa (), measurements of positively charged mesons are available from both the ALICE ALICE:2011aa (); Abelev:2012vra () and CDF Acosta:2003ax () collaborations. The latter sets of data offer new information on the charge separation of meson FFs that is not available from SIA where only the sum can be observed. It is also worth recalling that the Tevatron data from CDF Acosta:2003ax () are taken in rather than collisions and that LHCb has the unique capability to perform measurements at different asymmetric, forward rapidity intervals Aaij:2013mga (); Aaij:2015bpa (). Both sets of data will add unique information to our global analysis.

The factorized cross section for a given hadron and pseudorapidity may schematically be written as a convolution of appropriately combined PDFs, parton-to-hadron FFs, and partonic hard scattering cross sections:

 dσH1H2→hXdpTdη = 2pTS∑abcfH1a⊗fH2b⊗d^σcab⊗Dhc. (6)

Here, and denote the PDFs with flavor and in hadron and , respectively, and is the FF. The sum in (6) is over all contributing partonic cross sections , denoted as , which may be calculated as a perturbative series in , starting at which corresponds to the LO approximation. Hence, to perform a consistent NLO analysis of fragmentation functions, we include the corrections which have been computed analytically in ref:ppxsec-nlo (). As mentioned above, the factorized form given in Eq. (6) is again only valid up to power corrections that are suppressed by inverse powers of the hard scale, in this case . Throughout this work, we choose the factorization and renormalization scales for this process to be equal to the transverse momentum of the observed hadron, i.e., , but we will illustrate the residual dependence on in the phenomenological section below by varying by the conventional factor of two up and down.

One drawback of the single-inclusive high- production process is that the information on the dependence of the probed FFs is only accessible in integrated form through one of the convolution integrals in Eq. (6). The range of integration allowed by kinematical considerations for a given and of the observed hadron is rather broad and may reach well below the cut mentioned above. However, it has been shown in Ref. Sassot:2010bh () that one samples on average predominantly fairly large values of in Eq. (6), at mid rapidity and further increasing towards forward rapidities, and that values below are irrelevant for all practical purposes. Considering hadrons inside jets rather than single-inclusive hadron production allows one to sample more directly, as we shall discuss in some detail next.

### ii.3 D∗ meson in-jet production

The inclusive production of identified hadrons inside a fully reconstructed jet , where the hadron is part of the jet, has been studied for collisions in Refs. Arleo:2013tya (); Kaufmann:2015hma (); Kang:2016ehg (). In Arleo:2013tya (), the NLO cross section was obtained using a Monte-Carlo (MC) phase space integrator. Instead, in Refs. Kaufmann:2015hma (); Kang:2016ehg () analytical results were obtained using the approximation that the jet is sufficiently collimated. The NLO result of Kaufmann:2015hma () was derived within the standard pQCD framework, whereas Kang:2016ehg () employed methods within SCET for inclusive jet production Kang:2016mcy (); Dai:2016hz (), which allows for the additional resummation of single logarithms of the jet size parameter . At NLO accuracy, the analytical result of the cross section can be schematically written as . If the jet is sufficiently narrow, i.e., , power corrections of the order can be neglected. In studies for inclusive jet production ref:inclJets (), it was found that this “narrow jet approximation” is valid even for relatively large values of . For example, for the agreement between the thus obtained analytical results and the full MC result at NLO is better than 5%. This observation was also confirmed for the in-jet production of hadrons in Ref. Kaufmann:2015hma () by comparing to the full NLO MC calculation of Arleo:2013tya ().

In this work, we need the hadron-in-jet results for the anti- jet algorithm Cacciari:2008gp (). Currently, the only available data set for mesons within jets is provided by the ATLAS collaboration Aad:2011td () for which the anti- algorithm was used with a jet size parameter of . However, the results for cone Salam:2007xv () and Georgi:2014zwa () jets are also available in the literature Kaufmann:2015hma (); Kang:2016ehg (); Kang:2017mda ().

As it was discussed in detail in Ref. Kaufmann:2015hma (), the in-jet fragmentation provides a more direct access to the -dependence of FFs than data on single-inclusive hadron production. At LO accuracy, the cross section is directly proportional to the FFs probed at the momentum fraction , where

 zh≡pTpjetT (7)

and () denotes the transverse momentum of the hadron (jet). The cross section for the process may be written as

 dσpp→(jeth)XdpjetTdηjetdzh = 2pjetTS∑a,b,c∫1xminadxaxafa(xa,μ) (8) ×∫1xminbdxbxbfb(xb,μ)∫1zmincdzcz2cd^σcab(^s,^pT,^η,μ)vdvdw ×Ghc(zc,zh,μ,R),

where, again, we have set the renormalization and factorization scales to be equal and collectively denoted them by . For this process, we choose as our default choice of scale. The partonic cross sections are the same as they appear in the cross section for single-inclusive hadron production in Eq. (6). These hard functions depend on the jet partonic transverse momentum , the partonic rapidity and the partonic c.m.s. energy squared with the hadronic c.m.s. energy. The integration limits are customarily expressed in terms of the hadronic variables

 V≡1−pjetT√Se−ηjet,W≡(pjetT)2SV(1−V), (9)

 xmina = W,xminb=1−V1−VW/xa, zminc = 1−Vxb+VWxa. (10)

The function in Eq. (6) contains all the information on the production of the final-state jet and the identified hadron inside the jet and, hence, depends on the jet size parameter . To NLO accuracy, we can further decompose as

 Ghc(zc,zh,μ,R) = ∑ejc→e(zc,R,μ) (11) ×∑c′∫1zhdξξ~je→c′(ξ,R,μ)Dhc′(zhξ,μ).

The jet functions and describe the formation of the jet and the partonic fragmentation, respectively, and may be found in Ref. Kaufmann:2015hma (). Inserting Eq. (11) into Eq. (8), we may write the cross section as

 dσpp→(jeth)XdpjetTdη%jetdzh=∑e,c′Ee×[~je→c′⊗Dhc′](zh) (12)

where contains all the sums and integrals over the PDFs, the partonic cross sections and the jet functions and may be regarded as an “effective charge” weighting the different channels. The fragmentation functions appear in an actual convolution with the jet functions with respect to , multiplied by these effective charges. Eq. (12) illustrates the structural similarity of the in-jet fragmentation cross section and SIA, enabling access to the -dependence of the FFs. Due to the hadronic initial-state, the gluon fragmentation function already appears at LO accuracy, as it is the case for single-inclusive hadron production in collisions.

Typically, the hadron-in-jet production data are normalized to the inclusive jet cross section . Hence, the actual experimental observable is given by

 (13)

It was found in Kaufmann:2015hma (); Kang:2016mcy () that the cross section for inclusive jet production may be written in a similar form as the single-inclusive hadron production cross section in Eq. (6), with only the fragmentation functions replaced by perturbatively calculable jet functions , i.e.,

 dσH1H2→jetXdpjetTdη% jet=2pjetTS∑abcfH1a⊗fH2b⊗d^σcab⊗Jc. (14)

Thus, we may use the numerically efficient codes of Refs. ref:inclJets () to compute the hadron-in-jet cross section observable (13) in our global analysis of FFs.

## Iii Outline of the analysis

### iii.1 Parametrization

As we choose to work in the ZMVFNS for our global analysis of FFs, we closely follow the procedures for light hadron (pion and kaon) FFs as outlined in Refs. ref:dss (); Anderle:2015lqa (). However, due to the significantly smaller amount of data for production, we adopt a slightly less flexible, more economical functional form to parametrize the non-perturbative parton-to- FFs at some initial scale in the commonly adopted scheme:

 DD∗+i(z,μ20)=Nizαi(1−z)βiB[2+αi,βi+1]. (15)

We have tested that Eq. (15) nevertheless yields a very satisfactory description of the data, see also our results in Sec. IV below. The much simpler functional form with significantly less parameters also has the additional benefit of greatly facilitating the fitting procedure and the determination of uncertainties with the Hessian method.

We choose our initial scale to be equal to the charm quark mass . As we adopt the CT14 set of NLO PDFs and determination of the strong coupling  Dulat:2015mca () in all our calculations of hadronic cross sections and the scale evolution of FFs, we also use the heavy quark masses according to CT14, i.e., and . Furthermore, we assume that at the initial scale the FFs for all light quarks and antiquarks as well as for the anti-charm quark vanish, i.e.,

 DD∗+q(z,μ20)=0,       for q=u,¯u,d,¯d,s,¯s,¯c, (16)

which has no impact on the quality of the fit. In any case, none of these FFs can be reliably determined from the existing sets of data.

The bottom quark and antiquark FFs are included in the scale evolution above , and, as non-perturbative input, we only parametrize the total bottom-to- fragmentation function . In total this leaves us with 9 non-zero parameters in Eq. (15) for which is further reduced to 8 actual parameters to be determined in our global analysis since it turns out that is essentially unconstrained by data and has to be fixed.

As usual, the FFs for positively and negatively charged mesons are assumed to be related by charge conjugation, i.e.,

 DD∗−q(z,μ2)=DD∗+¯q(z,μ2) (17)

for quarks and

 DD∗−g(z,μ2)=DD∗+g(z,μ2) (18)

for the gluon. This will be used to compute cross sections for all data sets which observe only the sum of charges .

Finally, the parametrization in Eq. (15) is normalized to the respective Mellin moment by the denominator containing the Euler Beta function . Hence, the coefficients constitute the contribution of to the energy-momentum sum rule

 ∑h∫10dzzDhi(z,μ2)=1. (19)

### iii.2 Selection of data sets

Numerous experimental data exist for the three types of processes described in Sec. II. Identified mesons in collisions have been measured both by the ALEPH Barate:1999bg () and OPAL Akers:1994jc (); Ackerstaff:1997ki () collaborations at LEP at a c.m.s. energy of , the mass of the boson. Unfortunately, the results from the more recent OPAL analysis Ackerstaff:1997ki () are presented only in graphical form, and the corresponding numerical values are not anymore available ref:OPAL_pc (). Thus, we decide to use only the older set of OPAL data Akers:1994jc () with less statistics and somewhat larger uncertainties in our fit. Both collaborations also present bottom and charm flavor tagged data. Here, only the OPAL collaboration provides numerical values which we include in our global analysis.

Both data sets from ALEPH and OPAL are not corrected for the branching ratios of the decay channels used for the identification of the mesons. To obtain properly normalized cross sections, we divide the data by the branching ratios and which are given by Olive:2016xmw ()

 B1(D∗+→D0π+) = (67.7±0.5)%, B2(D0→K−π+) = (3.93±0.04)%. (20)

The uncertainties of the branching rations, and , are propagated into the systematic uncertainty of the SIA cross section data by adding them in quadrature, i.e.,

 Δdσsys= ⎷(dσΔB1B21B2)2+(dσΔB2B1B22)2+(ΔdσB1B2)2. (21)

Here, denotes the systematic error of the SIA data as provided by the ALEPH and OPAL experiments.

At lower c.m.s. energies, there are several measurements available around GeV ref:other-sia (). However, these data sets are rather old, and they consist of only a few data points that have very large uncertainties. Therefore, these sets do not add any relevant additional constraints to our global analysis, and, for simplicity, we choose to not include them.

Finally, some experiments have measured production just below the bottom threshold at around GeV. The most recent and precise data are from the BELLE collaboration Seuster:2005tr (). However, as stated on the HepData webpage, the “data for this record have been removed at the request of the authors due to an unrecoverable error in the measurement”; see ref:BELLEremoved (). Hence, we have to refrain from using this data set in our fit, which, potentially, could have been a very promising constraint from SIA in addition to the LEP data at . We note that the previous analysis of -meson FFs by the KKKS group Kneesch:2007ey (), to which we compare later on, includes the BELLE data as they were not yet withdrawn at the time when their fit was performed. Furthermore, CLEO Artuso:2004pj () and ARGUS Albrecht:1991ss () also provide data measured at similar c.m.s. energies as BELLE. However, both data sets are not corrected for initial-state radiation (ISR) effects. In addition, the ARGUS data points have very large uncertainties. The CLEO data have been included in the extractions of FFs in Refs. Cacciari:2005uk (); Kneesch:2007ey (). While Ref. Cacciari:2005uk () models the ISR effects based on data, the extraction of Kneesch:2007ey () includes ISR using certain approximations in the theory calculation of the cross section. However, both find noticeable tensions between the CLEO and ALEPH data. Since this may or may not be related to the treatment of ISR corrections, we choose not to include any of the low-energy SIA data in our analysis.

Data for inclusive production in hadronic collisions are available from the CDF collaboration at the Tevatron Acosta:2003ax () and from the ALICE ALICE:2011aa (); Abelev:2012vra (), ATLAS Aad:2015zix (), and LHCb Aaij:2013mga (); Aaij:2015bpa () collaborations at the LHC. We utilize all of these data sets in our global QCD analysis, as they provide valuable constraints on the gluon fragmentation function. As was mentioned in the Introduction and in Sec. II.3, an important new asset of our analysis are the in-jet fragmentation data for which ATLAS has presented results for identified mesons inside fully reconstructed jets Aad:2011td ().

To ensure the validity of the ZMVFNS approximation and the massless treatment of the mesons in the factorized formalism used to describe fragmentation processes, we have to impose certain cuts on the above mentioned data sets. For SIA, we only use data in the interval , i.e., , which is sufficient for the LEP data taken at . For all -spectra of mesons in hadronic collisions we select a very conservative cut of below which we exclude all data from the fit. Notice that this cut forces us to exclude the LHCb data sets from the 7 TeV run; we nevertheless show a comparison of our optimum fit to these data in Sec. IV. We are confident that our resulting set of FFs is not affected by our choice of since we find that lowering the cut down to does not lead to any significant changes in both the quality of the fit and the obtained optimum fit parameters in Eq. (15). This also implies that our results can be reliably extrapolated down to values of about , as we will also illustrate in some detail in Sec. IV.

### iii.3 Mellin Moment Technique

As mentioned above, we work entirely in complex Mellin moment space in order to solve the scale evolution equations of the FFs, to compute the relevant SIA and cross sections discussed in Sec. II, and to perform the actual fit and error analysis. The Mellin integral transform is well suited for these tasks as convolution integrals turn in ordinary products in Mellin space and the integro-differential evolution equations can be solved analytically. The resulting numerical codes for global QCD analyses are very efficient and fast.

In general, the pair of Mellin integral and inverse transforms of a function and are defined by

 f(N)=∫10dzzN−1f(z) (22)

and

 f(z)=12πi∫CNdNz−Nf(N), (23)

respectively, where denotes a suitable contour in the complex Mellin plane that guarantees fast convergence, see Refs. ref:pegasus (); Anderle:2016czy (); Anderle:2016kwa () for a comprehensive discussion of technical details and subtleties. In practice, one ends up having to compute only a limited number of moments along the contour in order to numerically solve the integral in Eq. (23).

Our analysis is set up in the following way: for each data point we use the analytical “truncated” solution of the evolution equations at NLO accuracy in Mellin space, see, e.g., Refs. ref:pegasus (); Anderle:2016czy (); Anderle:2016kwa (), to evolve the input FFs in (15) to the relevant scale. Next, the FFs are combined with appropriate space expressions for the hard scattering subprocesses before the inverse transform in Eq. (23) is performed numerically to evaluate the quality of the fit, see the next subsection. More specifically, in case of SIA, see Sec. II.1, this is achieved by taking the Mellin moments of Eq. (4) analytically, and all convolutions of FFs and coefficient functions turn schematically into

 D(z)⊗C(z)=12πi∫CNdNz−ND(N)C(N). (24)

Here, each coefficient function can be evaluated explicitly using the general definition in Eq. (22) and appropriate analytic continuations of harmonic sums to non-integer, complex values, see, for instance, Ref. Anderle:2016czy ().

For the more complicated expressions in scattering discussed in Secs. II.2 and II.3, one has to invoke an intermediate step as it is no longer possible or too cumbersome to perform the Mellin transform of the hard scattering cross sections analytically. Instead, we follow the steps outlined in Ref. Stratmann:2001pb (); ref:dss () and first express the FFs that appear, e.g., in Eq. (6) in terms of their respective Mellin inverse, see Eq. (23). After some reordering, the inclusive hadron production cross section can be recasted as follows

 dσH1H2→hXdpTdη = ∑c12πi∫CNdNDhc(N) ×2pTS∑abfH1a⊗fH2b⊗d^σcab⊗~Dhc,

where . The second line is independent of the FFs we are interested in, only needs to be evaluated once, and can be stored on a grid. In the end, for each data point one only has to perform the remaining contour integral in (III.3). This method is completely general and does not require any approximations such as -factors, and is also employed for the in-jet fragmentation in Sec. II.3.

### iii.4 Fitting and the Hessian Uncertainty Method

We obtain the optimum values for the eight free fit parameters in Eq. (15) by a standard minimization. We define the for the data sets included in the fit, each containing data points that pass the selection cuts specified in Sec. III.2, to be

 χ2=M∑i=1[(1−Ni)2ΔN2i+Mi∑j=1(NiTj−Ej)2ΔE2j], (26)

where is the experimental value for a given observable with uncertainty and is the corresponding theory calculation. Furthermore, we have introduced normalization shifts to account for this type of uncertainty whenever the normalization error is stated by the experiments. The optimum normalization shifts are computed analytically from the condition that they should minimize the . We note that we combine systematical and statistical uncertainties in quadrature in .

In order to estimate the uncertainties of the extracted FFs due to the experimental uncertainties and , we adopt the widely used iterative Hessian approach ref:hessian () to explore the range of possible variations of the obtained optimum parameters in the vicinity of the minimum of the function for a given tolerance . To this end, we provide 16 eigenvector sets for our FFs that correspond to the and directions of the eigenvectors of the diagonalized Hessian matrix. These sets greatly facilitate the propagation of hadronization uncertainties to any observable of interest. In fact, the uncertainty of an observable may be calculated straightforwardly as ref:hessian ()

 ΔO=12 ⎷8∑i=1[O+i−O−i]2, (27)

where denote the observable calculated with the plus or minus Hessian eigenvector set , respectively.

Finally, we note that choosing the tolerance is to some extent arbitrary in the presence of non-Gaussian or unaccounted uncertainties accompanying any global fit of PDFs or FFs. We have made sure that our Hessian sets, computed with , faithfully reflect the experimental uncertainties of the SIA data, as can be seen and will be discussed in the phenomenological section below. In this sense they correspond to uncertainties at the confidence level in the range that is constrained by data. Outside that range, the uncertainties are biased by the choice and flexibility of the selected functional form and assumptions made on the parameter space. In what follows, we will also briefly discuss additional, theoretical sources of uncertainties such as the choice and uncertainties of PDFs and from variations of the renormalization and factorization scales .

## Iv Results

### iv.1 Parton-to-D∗+ fragmentation functions

In this section, we present the results of our global determination of the parton-to--meson fragmentation functions and compare them to the previous fit of SIA data provided by Ref Kneesch:2007ey (), which will be labeled “KKKS08”. Note, that we use the public available numerical code from ref:KKKS08download (). In Tab. 1, we list the numerical values of the parameters of our optimal fit at NLO accuracy, see Eq. (15). As already mentioned, the parameter , which controls the behavior of the gluon FF, is basically unconstrained by data. For this reason, we decided to keep fixed. Note that other choices, like or 15 yield a total which differs by less than one unit, which is well within our tolerance . It is worth mentioning that in all fits with different values of , the parameter changes in such a way that the normalization remains essentially the same. As can be seen from the normalizations in Tab. 1, we find that the dominant contribution to mesons stems from valence charm quarks, , as is expected. The total bottom FF, contributes much less, and only a very small, though important, fraction of the gluon momentum is used to produce mesons. See the discussion of data below.

In Tab. 2 we list the data sets that pass the selection cuts on and as described in Sec. III.2 above and are thus included in our fit. We show the number of data points that are fitted for each set along with the obtained individual values. In addition, we present the analytically obtained optimum normalization shifts for each data set . They contribute to the total as specified in Eq. (26) and according to the quoted experimental normalization uncertainties ; an entry in Tab. 2 indicates that normalization uncertainties are not provided by the experiment. As can be seen, 96 data points from 3 different types of processes, SIA, single-inclusive hadron production, and in-jet fragmentation in collisions are included in our global QCD analysis of FFs for mesons, yielding a per degree of freedom of 1.17 for our best fit.

The so obtained FFs are shown for two representative scales and in Fig. 2 and 2, respectively, along with our uncertainty estimates (shaded bands) based on the Hessian method with ref:sets (), see Sec. III.4.

As the gluon and the unfavored light quark contributions turn out to be very small compared to the dominant charm-to- FF, we show them again for the sake of better legibility in the top right panel of Fig. 2. Notice that we just show the total FF as one example of the unfavored light quark and FFs, which are all the same as they are generated solely by QCD evolution from a vanishing input distribution, see Sec. III.1. This affects also the uncertainty estimates for these FFs which arise, again, just from evolution, i.e., mainly by propagating the uncertainties of the gluon FF. Hence, there is no direct access to the uncertainties of the unfavored light quark and FFs, such that they have to be taken with a grain of salt. Since none of the presently available data sets is sensitive to the unfavored FFs into a meson, in contrast to the also small gluon FF, one is forced to make some assumption about them. In any case, light quarks are expected to fragment mainly into light mesons such as pions and kaons, and their contribution to meson production should be small. Therefore, our choice of a vanishing input distribution for all unfavored FFs appears to be reasonable. We note that a similar assumption was made in the KKKS08 fit Kneesch:2007ey ().

It is instructive to compare the results of our FFs into mesons to those obtained in the KKKS08 fit Kneesch:2007ey () that is based only on SIA data and includes the by now obsolete and withdrawn BELLE data. The KKKS08 results for and the ratio to our FFs are shown as dashed lines in Figs. 2 and 2. As can be seen, one of the main differences is that our fit returns a significantly larger gluon contribution compared to KKKS08 at intermediate values of which might be related to the fact that also the gluon FF starts from a vanishing input in the KKKS08 fit. However, both the inclusive high- and, in particular, the in-jet fragmentation data, for the first time included in our global analysis, demand a non-zero gluon FF at our input scale in order to arrive at a satisfactory description of the data; see also the detailed comparisons to the inclusive and in-jet data below. One also notices, that the two valence charm FFs are somewhat shifted in with respect to each other and that also the height of the peak is different. This is most likely caused by the different sets of SIA data included in our and the KKKS08 analyses. Also, the KKKS08 fit does not include any uncertainty estimates.

Finally, in Fig. 2, for , we also show the bottom-to- FF which starts to evolve from a non-zero input above the threshold , see Sec. III.1. The total FF turns out to be quite similar to the one obtained in the KKKS08 analysis. This is to be expected as the bottom FF is largely constrained by the bottom-tagged data of the OPAL collaboration which are included in both fits.

### iv.2 Detailed comparison to data

In this section we compare theoretical calculations based on the results of our global QCD analysis with the available data. Throughout, we shall also show uncertainty bands obtained with the Hessian sets for as discussed in Sec. III.4. In addition, we perform all calculations with the FFs provided by Ref. Kneesch:2007ey (). Notice that Kneesch:2007ey () provides two sets of FFs which differ in the way they include finite charm quark and meson mass effects. Since we work in the ZMVFNS, we choose, as in Figs. 2 and 2 above, the corresponding KKKS08 set of FFs without quark mass effects in order to arrive at a meaningful comparison with our results. However, according to Kneesch:2007ey (), the KKKS analysis, some kinematic corrections due to the mass of the meson have been retained in all their fits beyond the standard theoretical framework based on factorization, which might be the source of some of the differences we observe at small .

We start with a study of the inclusive SIA data with identified mesons. The upper panel of Fig. 3 shows the LEP data at from the ALEPH Barate:1999bg () and OPAL Akers:1994jc () collaboration along with the theory calculations for the SIA multiplicities at NLO accuracy as defined in Eq. (1). The solid and dashed lines are obtained with our best fit, labeled as “this fit” throughout this section, and the FFs of KKKS08. The ratio of data over theory, our relative uncertainty estimates, and the ratio between KKKS08 and our best fit are given in the lower panel of Fig. 3. The hatched regions with and are excluded from our fit as discussed in Sec. III.2. The latter cut, which has no impact on the fit and in total only removes a single data point from our analysis, is imposed due to the presence of potentially large logarithms as , which cannot be properly accounted for in a fixed-order calculation.

As can be already anticipated from the individual values listed in Tab. 2, we find that our fit describes the inclusive SIA data very well, and our Hessian uncertainty estimates reflect the experimental uncertainties except for the data points with the lowest value of in each data set. Both sets of FFs describe the data equally well in the intermediate -region. Towards larger values of , the KKKS08 FFs overshoot the LEP SIA data significantly, which might be related to some tension with the CLEO and the by now withdrawn BELLE data, that are both included in their fit. In the small- region around our cut , the KKKS08 fit agrees slightly better with the data which might indicate some signs of a breakdown of the massless framework which we pursue in our analysis. In addition, we note, that the fixed-order evolution of FFs becomes more and more unstable towards smaller values of , see e.g. Anderle:2016czy (). Eventually, this can result in unphysical negative values for the FFs and the cross sections. The onset of this pathological behavior depends of the fit parameters and might, in part, also be responsible for the KKKS08 results to start to drop. At smaller values of than shown in Fig. 3, well below our cut , even our, still rising SIA multiplicity will start to drop and eventually reach unphysical, negative values.

In Fig. 4, we show the charm and bottom flavor-tagged data from OPAL Akers:1994jc () which are normalized to the total hadronic cross section. The bottom-tagged data are particularly instrumental in constraining the total bottom-to- FF in both our and the KKKS08 fit. As can be seen, both theoretical results describe the flavor-tagged data equally well, which have rather large uncertainties compared to the inclusive results shown in Fig. 3.

Following the order of processes as discussed in Sec. II, we next consider the single-inclusive, high- production of mesons in hadronic collisions. Since we are working in the ZMVFNS, we are especially interested in data where the observed meson has a transverse momentum much larger than the charm quark or the meson mass, i.e., . As discussed in Sec. III.2, we employ a rather stringent cut of in our global analysis. However, we will demonstrate that the so obtained FFs work unexpectedly well in describing single-inclusive meson cross sections down to much smaller values of around .

In this respect, the most relevant data set is the one presented by the ATLAS collaboration Aad:2015zix (), shown in Fig. 5, which covers the range at a c.m.s. energy of integrated over the mid rapidity range . In the upper panel, a comparison of the ATLAS data with calculations at NLO accuracy is presented based on our best fit and KKKS08 FFs and using Eq. (6); again, data in the hatched area, i.e., below our cut , data are not included in our global analysis. The middle panel gives the ratios of the KKKS08 prediction and the ATLAS data with respect to our NLO calculation. In addition, it illustrates the uncertainty estimates (shaded bands) obtained from our Hessian sets of FFs.

Both sets of FFs provide a satisfactory description of the data at NLO accuracy in the ZMVFNS even well below . About of the mesons at originate from gluon fragmentation which drops down to approximately at the highest measured by ATLAS. In view of the sizable differences between our and the KKKS08 gluon FF illustrated in Figs. 2 and 2, the similarity of the NLO cross sections is a remarkable result and indicates that the inclusive spectra only constrain certain -moments of the gluon FF rather than its detailed shape. The differences between the two sets of FFs in Figs. 2 and 2, in particular, the gluon FF, will be much more pronounced when we turn to the in-jet fragmentation data below.

In the lower panel of Fig. 5, we show other important sources of theoretical uncertainties associated with a pQCD calculation of cross sections based on, e.g., Eq. (6). The outer shaded bands illustrates the ambiguities due to simultaneous variations of the factorization and renormalization scales in the range . As can be seen, in the range used in our fit, this results in a roughly constant relative uncertainty of about . The theoretical error related to PDF uncertainties are estimated with the Hessian sets provided by the CT14 collaboration Dulat:2015mca (). They turn out to be smaller than QCD scale uncertainties and are at a level of about as can be inferred from the inner shaded bands in the lower panel of Fig. 5. To be compatible with our estimates of the one-sigma uncertainties of the FFs we follow Ref. Watt:2011kp () and rescale the available CT14 Hessian sets from the to the confidence level by a applying constant factor .

A large amount of data points for inclusive -meson production have been presented by the LHCb collaboration. They measured the single-inclusive production cross section at forward rapidities for two different c.m.s. energies,  TeV Aaij:2013mga () and  Aaij:2015bpa (). For each c.m.s. energy, the data are presented in five bins of rapidity in the range from up to . Compared to the mid rapidity data shown in Fig. 5, the LHCb data are limited to smaller values of . Nevertheless, several data points from the run Aaij:2015bpa () are above our cut except for the most forward rapidity bin but, unfortunately, none of the data points taken at  Aaij:2013mga () passes the cut. Both sets of data are shown in Figs. 6 and 7 and compared to the results of NLO calculations based on our and the KKKS08 set of FFs. We note that the more forward the rapidity interval, the more important is the role of gluon fragmentation in producing the observed mesons, a feature that has already been observed for the production of lighter hadrons at the LHC Sassot:2010bh (). For instance, at around of the mesons at originate from gluons. Since forward data also sample on average larger values of Sassot:2010bh (), the LHCb data nicely complement the mid rapidity data by ATLAS.

As for the ATLAS data, both sets of FFs also give an equally good description of the LHCb data shown in Fig. 6 for , as can be also inferred from Tab. 2, and they continue to follow the data well below our cut, down to about . Also the data taken at , that are not included in our fit, are well described down to except for the most forward bin . The KKKS08 FFs follow the trend of the data even further down to the lowest values shown in Figs. 6 and 7; for the sake of applicability of pQCD, we refrain from showing comparisons to the LHCb data below . This feature of the KKKS08 fit, which is unexpected in a ZMVFNS approach, might be due to the inclusion of finite hadron mass corrections in their fit of SIA data, that are, however, beyond the factorized framework outlined in Sec. II and adopted by us. It is also interesting to notice that there are some indications for a mild tension between the ATLAS and the LHCb data in our global fit. The ATLAS data alone would prefer a somewhat larger gluon-to- meson FF as can be inferred from the middle panel of Fig. 5. This would yield a significantly better fit of the ATLAS data in terms of even when the in-jet fragmentation data, which we shall discuss next, are included in the fit. The latest, revised version of the LHCb data Aaij:2015bpa () does not tolerate, however, such an increased gluon FF in our global analysis.

We refrain from showing comparisons of our theoretical results with the ALICE and CDF data on single-inclusive, high- meson production. As can be seen from Tab. 2, the few data points which pass our cut on are very well reproduced by our fit. Again, adopting the KKKS08 set of FFs leads to a similar description of these data, assuming .

Finally, we turn to data on in-jet production, which, in this paper, are considered for the first time in a global QCD analysis of FFs and, hence, represent the centerpiece of our phenomenological studies. The relevant QCD formalism to compute in-jet production in the standard factorized framework at NLO accuracy was sketched in Sec. II.3. The main and novel asset of this process, as compared to single-inclusive hadron production in collisions, is the fact that in-jet data probe the parton-to-hadron FFs locally in the momentum fraction in the LO approximation. Therefore, one anticipates a much improved sensitivity to the, in particular, -dependence of the gluon FF also beyond LO accuracy than from single-inclusive probes. In the latter case, we have just found that two rather different gluon FFs, ours and the one from the KKKS08 fit, can result both in a good description of the existing data, cf. Figs. 5 and 6 above.

Specifically, for the in-jet production of mesons, it was found in Ref. Chien:2015ctp () that the cross section computed with the KKKS08 set of FFs falls significantly short of the corresponding yields observed by ATLAS Aad:2011td (). The authors of Ref. Chien:2015ctp () observed that by increasing the KKKS08 gluon FFs ad hoc by a -independent factor of 2 would help to better describe the ATLAS data. However, such a modified gluon FF would then significantly overshoot the single-inclusive data for -mesons. Clearly, to address this issue reliably and in detail, a simultaneous global QCD analysis of all relevant probes comprising SIA, and single-inclusive and in-jet production in proton-proton collisions is absolutely essential.

From Fig. 8 and Tab. 2 one can gather that our global fit yields an excellent description of the in-jet data by ATLAS in all five bins of the jet’s transverse momentum without compromising the comparison to SIA or single-inclusive data. A corresponding calculation with the KKKS08 set of FFs falls short of the data for momentum fractions of the meson, as was already observed in Ref. Chien:2015ctp (). In fact, the -dependence of the NLO calculations with the two different sets of FFs very closely follow the corresponding dependence of the gluon-to- FF illustrated for two different scales in Figs. 2 and 2. The main difference between our analysis and the KKKS08 extraction of FFs is that we allow for a non-zero gluon FF at our initial scale which appears to be necessary in order to achieve a good global fit of all data; recall that the KKKS08 analysis was based only on SIA data where some assumption about the gluon FF has to be made. The quark FFs, in particular, the charm FF, adjust accordingly in the fit but play only a very minor role in computations of cross sections in the range currently covered by experiment. Finally, we note that theoretical uncertainties due to the choice of scale and from ambiguities in the adopted set of PDFs are of similar size as we have estimated for the single-inclusive data; cf. the lower panel of Fig. 5 and Ref. Kaufmann:2015hma ().

Our case study of clearly reveals how powerful in-jet data can be in further constraining FFs. Based on the framework developed and applied in this paper, in-jet data can be straightforwardly included in any future global fit of FFs once such data become available.

## V Conclusions and Outlook

We have presented the first global QCD analysis of fragmentation functions that makes use of in-jet data besides the usual sets of experimental results stemming from single-inclusive hadron production in electron-positron annihilation and proton-(anti)proton collisions. The necessary technical framework to incorporate in-jet fragmentation data consistently into a global fit at next-to-leading order accuracy was outlined in detail, and an implementation within the Mellin moment technique was given and henceforth adopted in all our phenomenological studies.

As a case study, we have analyzed available data for charged mesons in terms of parton-to- meson fragmentation functions. An excellent global description of all the different processes included in the fit was achieved. In particular, the in-jet fragmentation data have been shown to be of great importance in pinning down the otherwise largely unconstrained momentum fraction dependence of the gluon fragmentation function. Compared to the only other previously available set of fragmentation function, that was based solely on electron-positron annihilation data, we obtain a rather different momentum dependence for the hadronization of gluons in order to describe the in-jet data.

In addition to our optimum fit, we have also, for the first time, estimated the uncertainties of charged meson fragmentation functions. To this end, we have applied the Hessian method. The obtained Hessian sets provide a straightforward way to propagate our estimated uncertainties to any other process of interest. Apart from the experimental uncertainties that are incorporated in the Hessian sets, we have illustrated the importance of other, theoretical sources of ambiguities comprising the actual choice of renormalization and factorization scales and corresponding uncertainties of parton distribution functions which are needed in calculations of any hadronic collision process.

For the time being, we have adopted the ZMVFNS throughout our global analysis, i.e., we have imposed rather stringent cuts on the minimum transverse momentum of the mesons for data to be included in our fit. We have demonstrated, however, that our fit gives a reasonable description of single-inclusive data from the LHC both at mid and forward rapidities even down to significantly smaller values of transverse momentum of about .

We believe that in-jet data will prove very valuable in the future in any upcoming analysis of fragmentation functions, in particular, in further constraining the detailed momentum dependence of the hadronization of gluons. The framework developed and applied in this paper can be straightforwardly generalized to incorporate in-jet data in any future global fit of FFs once such data become available. We plan to extend our phenomenological studies to charged and neutral mesons in the near future. We note that the theoretical framework for the in-jet production of hadrons at next-to-leading order accuracy has been recently extended to include also photons Kaufmann:2016nux (). The fragmentation into photons is so far only rather poorly understood and constrained by data. Again, we expect any upcoming in-jet data to be very valuable in a new extraction of photon fragmentation functions. Finally, we plan to study in detail the impact of the resummation of logarithms of the jet size parameter . By making use of the results for the in-jet fragmentation of hadrons derived within the SCET formalism, it is possible to extract fragmentation functions at a combined accuracy of NLO+NLL.

## Acknowledgments

We are grateful to Y.-T. Chien, Z.-B. Kang, E. Mereghetti, N. Sato, and W. Vogelsang for helpful discussions and comments. D.P.A. was supported by the Lancaster-Manchester-Sheffield Consortium for Fundamental Physics under STFC Grant No. ST/L000520/1. T.K. was supported by the Bundesministerium für Bildung und Forschung (BMBF) under grant no. 05P15VTCA1. F.R. is supported by the LDRD Program of Lawrence Berkeley National Laboratory, the U.S. Department of Energy, Office of Science and Office of Nuclear Physics under contract number DE-AC02-05CH11231. I.V. is supported by the US Department of Energy, Office of Science under Contract No. DE-AC52-06NA25396 and the DOE Early Career Program.

### References

1. See, e.g., J. C. Collins, D. E. Soper, and G. F. Sterman, Adv. Ser. Direct. High Energy Phys. 5, 1 (1988).
2. A. Metz and A. Vossen, Prog. Part. Nucl. Phys. 91, 136 (2016).
3. R. Gauld, J. Rojo, L. Rottoli, and J. Talbert, JHEP 1511, 009 (2015).
4. M. V. Garzelli, S. Moch, and G. Sigl, JHEP 1510, 115 (2015); A. Bhattacharya, R. Enberg, M. H. Reno, I. Sarcevic, and A. Stasto, JHEP 1506, 110 (2015); R. Gauld, J. Rojo, L. Rottoli, S. Sarkar, and J. Talbert, JHEP 1602, 130 (2016); M. Benzke, M. V. Garzelli, B. Kniehl, G. Kramer, S. Moch, and G. Sigl, arXiv:1705.10386 [hep-ph].
5. Z. B. Kang, F. Ringer, and I. Vitev, JHEP 1703, 146 (2017); for a recent review, see, e.g., A. Andronic et al., Eur. Phys. J. C 76, 107 (2016).
6. W. Beenakker, H. Kuijf, W. L. van Neerven, and J. Smith, Phys. Rev. D 40, 54 (1989); P. Nason, S. Dawson, and R. K. Ellis, Nucl. Phys. B 327, 49 (1989); W. Beenakker, W. L. van Neerven, R. Meng, G. A. Schuler, and J. Smith, Nucl. Phys. B 351, 507 (1991). M. L. Mangano, P. Nason, and G. Ridolfi, Nucl. Phys. B 373, 295 (1992).
7. M. Cacciari and M. Greco, Nucl. Phys. B 421, 530 (1994); M. Cacciari, M. Greco, and P. Nason, JHEP 9805, 007 (1998); M. Cacciari, S. Frixione, and P. Nason, JHEP 0103, 006 (2001)
8. M. Cacciari, P. Nason, and C. Oleari, JHEP 0604, 006 (2006).
9. J. Binnewies, B. A. Kniehl, and G. Kramer, Phys. Rev. D 58, 014014 (1998); B. A. Kniehl and G. Kramer, Phys. Rev. D 74, 037502 (2006); B. A. Kniehl, G. Kramer, I. Schienbein, and H. Spiesberger, Phys. Rev. D 71, 014018 (2005); Eur. Phys. J. C 41, 199 (2005); Eur. Phys. J. C 72, 2082 (2012); Eur. Phys. J. C 75, no. 3, 140 (2015).
10. T. Kneesch, B. A. Kniehl, G. Kramer, and I. Schienbein, Nucl. Phys. B 799, 34 (2008).
11. M. Fickinger, S. Fleming, C. Kim, and E. Mereghetti, JHEP 1611, 095 (2016).
12. S. Kretzer, Phys. Rev. D 62, 054001 (2000); S. Albino, B. A. Kniehl, and G. Kramer, Nucl. Phys. B 725, 181 (2005); M. Hirai, S. Kumano, T.-H. Nagai, and K. Sudoh, Phys. Rev. D 75, 094009 (2007); N. Sato, J. J. Ethier, W. Melnitchouk, M. Hirai, S. Kumano, and A. Accardi, Phys. Rev. D 94, 114004 (2016).
13. D. de Florian, R. Sassot, and M. Stratmann, Phys. Rev. D 75, 114010 (2007); Phys. Rev. D 76, 074033 (2007); M. Epele, R. Llubaroff, R. Sassot, and M. Stratmann, Phys. Rev. D 86, 074028 (2012); D. de Florian, M. Epele, R. J. Hernandez-Pinto, R. Sassot and M. Stratmann, Phys. Rev. D 91, 014035 (2015); Phys. Rev. D 95, 094019 (2017).
14. M. Stratmann and W. Vogelsang, Phys. Rev. D 64, 114007 (2001).
15. D. P. Anderle, F. Ringer, and M. Stratmann, Phys. Rev. D 92, 114017 (2015).
16. V. Bertone, S. Carrazza, N. P. Hartland, E. R. Nocera and J. Rojo, arXiv:1706.07049 [hep-ph].
17. D. P. Anderle, F. Ringer, and W. Vogelsang, Phys. Rev. D 87, 034014 (2013).
18. D. P. Anderle, T. Kaufmann, M. Stratmann, and F. Ringer, Phys. Rev. D 95, 054003 (2017).
19. D. Anderle, D. de Florian, and Y. Rotstein Habarnau, Phys. Rev. D 95, 034027 (2017).
20. J. Pumplin, D. R. Stump, and W. K. Tung, Phys. Rev. D 65, 014011 (2001); J. Pumplin, D. Stump, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai, and W. K. Tung, Phys. Rev. D 65, 014013 (2001)
21. M. Procura and I. W. Stewart, Phys. Rev. D 81, 074009 (2010), Erratum: [Phys. Rev. D 83, 039902 (2011)]; A. Jain, M. Procura, and W. J. Waalewijn, JHEP 1105, 035 (2011); JHEP 1204, 132 (2012); C. W. Bauer and E. Mereghetti, JHEP 1404, 051 (2014).
22. F. Arleo, M. Fontannaz, J. P. Guillet, and C. L. Nguyen, JHEP 1404, 147 (2014).
23. T. Kaufmann, A. Mukherjee, and W. Vogelsang, Phys. Rev. D 92, 054015 (2015).
24. Z. B. Kang, F. Ringer, and I. Vitev, JHEP 1611, 155 (2016).
25. Y. T. Chien, Z. B. Kang, F. Ringer, I. Vitev, and H. Xing, JHEP 1605, 125 (2016).
26. G. Aad et al. [ATLAS Collaboration], Phys. Rev. D 85, 052005 (2012).
27. R. Bain, L. Dai, A. Hornig, A. K. Leibovich, Y. Makris, and T. Mehen, JHEP 1606, 121 (2016).
28. G. Aad et al. [ATLAS Collaboration], Phys. Lett. B 739, 320 (2014); B. Abelev et al. [ALICE Collaboration], JHEP 1403, 013 (2014); S. Chatrchyan et al. [CMS Collaboration], Phys. Rev. C 90, 024908 (2014); The ATLAS Collaboration, ATLAS-CONF-2017-004.
29. R. Aaij et al. [LHCb Collaboration], Phys. Rev. Lett. 118, 192001 (2017).
30. H. Abramowicz et al. [H1 and ZEUS Collaborations], JHEP 1509, 149 (2015).
31. P. Nason and B. R. Webber, Nucl. Phys. B 421, 473 (1994); ibid. B 480, 755 (1996) (E).
32. D. de Florian, M. Stratmann, and W. Vogelsang, Phys. Rev. D 57, 5811 (1998).
33. P. J. Rijken and W. L. van Neerven, Phys. Lett. B 386, 422 (1996); Nucl. Phys. B 487, 233 (1997); A. Mitov and S. O. Moch, Nucl. Phys. B 751, 18 (2006).
34. G. Aad et al. [ATLAS Collaboration], Nucl. Phys. B 907, 717 (2016).
35. R. Aaij et al. [LHCb Collaboration], Nucl. Phys. B 871, 1 (2013).
36. R. Aaij et al. [LHCb Collaboration], JHEP 1603, 159 (2016), Erratum: [JHEP 1609, 013 (2016)], Erratum: [JHEP 1705, 074 (2017)].
37. B. Abelev et al. [ALICE Collaboration], JHEP 1201, 128 (2012).
38. B. Abelev et al. [ALICE Collaboration], JHEP 1207, 191 (2012).
39. D. Acosta et al. [CDF Collaboration], Phys. Rev. Lett. 91, 241804 (2003).
40. F. Aversa, P. Chiappetta, M. Greco, and J. P. Guillet, Nucl. Phys. B 327, 105 (1989); B. Jager, A. Schafer, M. Stratmann, and W. Vogelsang, Phys. Rev. D 67, 054005 (2003).
41. R. Sassot, P. Zurita, and M. Stratmann, Phys. Rev. D 82, 074011 (2010).
42. Z. B. Kang, F. Ringer, and I. Vitev, JHEP 1610, 125 (2016).
43. L. Dai, C. Kim, and A. K. Leibovich, Phys. Rev. D 94, 114023 (2016).
44. F. Aversa, M. Greco, P. Chiappetta, and J. P. Guillet, Z. Phys. C 46, 253 (1990); B. Jager, M. Stratmann, and W. Vogelsang, Phys. Rev. D 70, 034010 (2004); A. Mukherjee and W. Vogelsang, Phys. Rev. D 86, 094009 (2012); T. Kaufmann, A. Mukherjee, and W. Vogelsang, Phys. Rev. D 91, no. 3, 034001 (2015).
45. M. Cacciari, G. P. Salam, and G. Soyez, JHEP 0804, 063 (2008).
46. G. P. Salam and G. Soyez, JHEP 0705, 086 (2007).
47. H. Georgi, arXiv:1408.1161 [hep-ph]; Y. Bai, Z. Han, and R. Lu, JHEP 1503, 102 (2015).
48. Z. B. Kang, F. Ringer, and W. J. Waalewijn, arXiv:1705.05375 [hep-ph].
49. S. Dulat et al., Phys. Rev. D 93, 033006 (2016).
50. R. Barate et al. [ALEPH Collaboration], Eur. Phys. J. C 16, 597 (2000).
51. R. Akers et al. [OPAL Collaboration], Z. Phys. C 67, 27 (1995).
52. K. Ackerstaff et al. [OPAL Collaboration], Eur. Phys. J. C 1, 439 (1998).
53. S. Söldner-Rembold, private communication.
54. C. Patrignani et al. [Particle Data Group], Chin. Phys. C 40, 100001 (2016).
55. J. M. Yelton et al., Phys. Rev. Lett. 49, 430 (1982); P. S. Baringer et al., Phys. Lett. B 206, 551 (1988); W. Bartel et al. [JADE Collaboration], Phys. Lett. 146B, 121 (1984); H. Aihara et al. [TPC/Two Gamma Collaboration], Phys. Rev. D 34, 1945 (1986); W. Braunschweig et al. [TASSO Collaboration], Z. Phys. C 44, 365 (1989).
56. R. Seuster et al. [Belle Collaboration], Phys. Rev. D 73, 032002 (2006).
57. Concerning the withdrawal of the BELLE data, see: http://durpdg.dur.ac.uk/view/ins686014.
58. M. Artuso et al. [CLEO Collaboration], Phys. Rev. D 70, 112001 (2004).
59. H. Albrecht et al. [ARGUS Collaboration], Z. Phys. C 52, 353 (1991).
60. A. Vogt, Comput. Phys. Commun. 170, 65 (2005).
62. A Fortran code of our optimum fit of , and FFs and the corresponding Hessian uncertainty sets are available upon request from the authors.
63. G. Watt, JHEP 1109, 069 (2011).
64. T. Kaufmann, A. Mukherjee, and W. Vogelsang, Phys. Rev. D 93, 114021 (2016).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters