Impact of the NOMAD data on the strangeness in the proton

(1)

T

HESIS

MS

C

P

ROJECT

T

HEORETICAL

P

HYSICS

Impact of the NOMAD data on the

strangeness in the proton

BY

FERRAN

FAURA

I

GLESIAS

Student nr. 11045914 Size 60 EC.

Conducted between September 2 2019 and July 17 2020. University of Amsterdam & Vrije Universiteit Amsterdam

NATIONAL

I

NSTITUTE FOR

S

UBATOMIC

P

HYSICS

Supervisor dr. J.ROJO

Examiner prof.dr. P.J.MULDERS

(2)

Abstract

Parton distribution functions (PDFs) are a vital part of high-energy physics phenomenology. Among all light quarks, the strange sea PDFs are the least known. Such a poor knowledge is aggravated by the fact that different determinations of the strange sea PDFs can be obtained if different scattering processes sensitive to it are included in a global PDF fit, namely neu-trino induced deep inelastic scattering (DIS) or gauge boson production in proton–proton collisions. To clarify this state of affairs, in this thesis we analyze the DIS measurements performed by the NOMAD experiment in the framework of the NNPDF3.1 PDF set for the first time. We find that they lead to a reduction of the strange sea PDF uncertainties and that they provide evidence for a suppression of the strange sea in comparison to the other light sea quarks.

(3)

1. Introduction

Improving our knowledge of the internal structure of the proton is an important part of the phenomenology programs related to high-energy particle physics. In any high-energy scattering experiment involving a hadron one measures a hadronic observable. We want to describe these observables in terms of fundamental partonic properties of the measured hadron. This description is realized by using parton distribution functions. These encode the momentum distributions of the quarks and gluons that make up the measured hadron. To this day there are no reliable theoretical approaches to compute PDFs from first princi-ples. This means that PDFs are determined from experimental data by solving a regression problem. There is a large variety of experimental data gathered by different collaborations, as different types of observables can be measured.

In this study we are interested in the strangeness in the proton, or i.e. the strange sea PDFs in the proton. The strange and anti-strange quark PDFs are the least known of all light quarks, due to the limited experimental data available that is sensitive to the strange content of the proton. The strangeness is commonly probed by the neutrino DIS process, νs→ cX, which in turn can semi-leptonically decay in the so-called dimuon final state as measured by the CCFR [1], NuTeV [2] and NOMAD [3] collaborations.

The reason for our interest in the strange content of the proton is due to opposing results of the strange sea suppression factor, which is given by [4]

Ks= s

(x, Q2) + ¯s(x, Q2)

¯

u(x, Q2_{) +}_d¯₍_{x, Q}2₎. (1.1)

The value of this quantity is found to be 0.55±0.21, 0.57±0.17, 0.60±0.13 and 0.63±0.03 for CT14 [5], MMHT2014 [6], NNPDF3.1 [7] and ABMP16 [8] computed at x = 0.023 and Q2 = 1.9 GeV. These values show a similar order of magnitude of Ks ' 0.5, which means

the strange content is suppressed relative to the non-strange light quark content. There is, however, an important opposing result measured to be 1.13+₋0.08_0.13 coming from an analysis done by the ATLAS collaboration based on the collider W,Z differential cross section data [9]. Due to the importance of strangeness for example in the measurement of the mass of the W boson, we make an attempt in this thesis to resolve this issue.

This project will focus on investigating the inclusion of the NOMAD dimuon data into the NNPDF3.1 set. The NOMAD dimuon data is based on a significantly larger sample size

(6)

CHAPTER 1. INTRODUCTION

compared to other dimuon experiments, therefore increasing its precision and potentially delivering a more precise determination of the strange sea. In general when fitting a PDF to a set of experimental data one chooses a parametrization, a method to represent the un-certainty of the data into the unun-certainty of the PDFs and finally a method to optimize the parameters to the data. NNPDF uses a neural network for the parametrization and a genetic algorithm to optimize the parameters. This neural network is provided with a set of Monte Carlo sampled data, which is then used to represent the uncertainty of the PDF. This process is in general complicated and can be severely time consuming. A useful method to perform a fit and avoid this loss of time is by means of reweighting. This method provides similar results as fitting new data, while consuming less time at the cost of losing information. It is, however, possible to quantify this loss of information to determine whether it is statistically significant.

To start with, in chapter 2 we will discuss the basic concepts of DIS and define the DIS kinematics. A simple derivation of a DIS cross section is given. In the following section there is an explanation of the Parton Model, asymptotic freedom, the QCD improved Parton Model that includes the momentum evolution of the PDFs and a definition of the structure functions.

In chapter 3 we give a brief explanation on the method of parametrization of the PDFs used by NNPDF, which is by means of a neural network. Furthermore we explain two sepa-rate methods for error propagation and finally we give a short summary on the strangeness-sensitive datasets included in the NNPDF3.1 set.

In chapter 4 we explain in detail the reweighting and unweighting procedure based on Refs. [10] and [11]. We discuss the quantification of the loss of information as a consequence of reweighting and a method to determine the consistency of the data.

In chapter 5 we explain the NOMAD observable based on Ref. [3]. First we explain the general process charm dimuon production which is measured by NOMAD. Then we give a general explanation of the observable. As the observable is a ratio between the charm dimuon cross section with the inclusive DIS neutrino-nucleon cross section, the remaining part of the chapter is divided in two parts. The numerator, which is the charm dimuon cross section, and the denominator, which is the inclusive DIS neutrino-nucleon cross section.

In chapter 6 we give the method used to reproduce the observable. The calculations are performed in a standalone PYTHON program which is benchmarked against results of Ref. [12]. Then we discuss the precise expressions implemented in the code, including the kine-matic limits which determine the limits of integration. Finally, we briefly summarize the methods used to evaluate the structure functions up to NNLO.

(7)

observ-CHAPTER 1. INTRODUCTION

able by means of both the global and the collider set. The global set is a PDF set which is fitted to the full spectrum of available types of data. This is in contrast to the collider set, which contains solely data gathered by collider experiments. In this chapter we show the dif-ference in strangeness determined by both these types of sets. Furthermore we compute the NOMAD observable at different values for Qmin to determine the contribution to the cross

section prediction of the kinematic region which is unreachable by NNPDF3.1. For compar-ison we compute the NOMAD observable using the CT18 and ABM12 PDF sets, which use an alternative method to compute the uncertainties. In this chapter we also show the results of the reweighting procedure, including the impact of reweighting on the PDFs. Finally we build an unweighted set, which is used to compute the strange sea suppression factor.

(8)

2. Deep Inelastic Scattering

A method to probe the internal structure of the proton is to bombard the proton with a beam of high energy particles. One of the breakthrough experiments showing that the proton was not a point-like particle was a DIS experiment done by SLAC in the 1960s. They used a 20 GeV electron beam incident on static proton targets, which was sufficiently energetic to reveal the internal structure of the proton. They showed that electrons were scattering off quasi-free point-like particles inside the proton which were later identified as quarks. DIS experiments have since been very important as it is one of the processes that provide detailed information on the PDFs of the internal constituents of the nucleons. In this section we describe the basics of DIS relevant to this study. First, we define the relevant kinematics related to DIS, then we derive the differential cross section that is studied in this project. Furthermore, we explain the Parton Model, the QCD improved Parton Model and finally give the definitions of the structure functions.

2.1 DIS Kinematics

We can describe a DIS process in a general way by considering inclusive lepton–nucleon scattering

ν(k) +N(p) → `(k0) +X(pX), (2.1)

where X is the unknown hadronic final state, ν, N and ` are the incoming neutrino and nucleon and the outgoing lepton with k, p and pX being the corresponding four-momenta.

In general we can distinguish two types of neutrino induced scattering: neutral current (NC) and charged current (CC) scattering. The first one is mediated by the exchange of a Z boson, while the latter is mediated by the charged W±boson. NOMAD considers incoming muon neutrinos νµ and a W+ exchange boson, therefore we will restrict ourselves to this type of

(9)

CHAPTER 2. DEEP INELASTIC SCATTERING

Figure 2.1: A basic DIS Feynman diagram with incoming neutrino with four-momentum k, outgoing lepton with momentum k0, the momentum transfer q = k−k0 carried by the exchange electroweak boson, the static proton momentum p and the undetected hadronic state X. [13]

First we start with the general process given in Eq. (2.1), which is depicted by the di-agram in figure 2.1. In this didi-agram, k and k0 are the incoming and outgoing lepton four-momenta, respectively. The proton in the initial state has four-momentum p, while the inclusive hadronic final state has four-momentum X. The exchange electroweak boson’s four-momentum is denoted by q and is equal to the difference of the incoming and outgoing lepton’s momenta, thus q = k−k0. This quantity is usually referred to as the momentum transfer. From these four momenta we can build a set of Lorentz invariant quantities which are commonly used to describe DIS kinematics:

Q2 = −q2 = −(k−k0)2 Momentum transfer, (2.2)

M2 = p2 Proton mass squared, (2.3)

s = (p+k)2 Centre-of-mass energy squared, (2.4)

ν= p·q

M =E−E

0

Energy transfer in the laboratory frame, (2.5) x = Q

2

2p·q = Q2

2Mν Bjorken scaling variable, (2.6) y= p·q

p·k =

ν

E Inelasticity; fractional energy transfer, (2.7) W2 = (p+q)2 = M2−Q2+2Mν Invariant mass of X squared, (2.8) where M is the proton mass and E(E0)is the energy of the incoming(outgoing) lepton. We shall see later that the Bjorken scaling variable x corresponds to the fraction of proton longi-tudinal momentum carried by the struck quark in the reference frame in which the proton has infinite momentum. As these quantities are Lorentz invariant, we can choose a conve-nient frame to calculate them. The four momenta are then defined in the laboratory frame

(10)

as follows:

k = (E,~k) (2.9)

k0 = (E0,~k0) (2.10)

p = (M,~0) (2.11)

The process is called deep and inelastic due to the limits Q2 M2 _{and W}2 _M2_{. This}

means that since we are in the deeply inelastic region we can neglect the masses of the initial and final state leptons:

k2=k02≈0, (2.12) from which it follows that

E≈ |k|, E0 ≈ |k0|. (2.13)

2.2 Differential cross section

The differential cross section dσ per unit solid angle dΩ is given by dσ dΩ = 1 64π2_s |p0| |p||M| 2_, _(2.14)

with s being the square of the centre-of-mass energy, M the scattering amplitude and|p|

and|p0|the absolute values of the three-momenta of the initial and final state respectively. If all four particles in the interaction have equal mass, this expression simplifies to

dσ dΩ =

|M|2

64π2_s. (2.15)

We want to write the expression for an unpolarized inclusive neutrino–nucleon(νN → `X)

DIS process. For this we use (2.15), which is the general expression for the differential cross section.

To calculate the cross section we have to construct the square of the scattering amplitude

|M|2_{. In principle one would think that we can build} _|M|2 _{by using the QED Feynman}

rules related to the diagram of figure 2.1. However this is not possible, because there are no Feynman rules for the hadronic vertex. This means that we have to define a tensor function, Wµν, that represents the unknown hadronic physics. However, first we write the scattering

amplitude and its Hermitian conjugate as iM(νN→ `X) =i√GF 2 ¯ `(k0)γµ(1−γ 5₎ 2 ν(k) hX|Jµ(q) |Pi, (2.16) −iM†₍ νN→ `X) = −i√GF 2¯ν(k)γ ν(1−γ 5₎ 2 `(k 0_{) h}_P_|_J† ν(q) |Xi, (2.17)

where GF is the Fermi coupling that gives the strength of the weak interactions [14], Jµ(q)

is the electroweak current at the hadronic vertex and γµ₍₁₋

γ5)is the helicity contribution

(11)

CHAPTER 2. DEEP INELASTIC SCATTERING |M(νN→ `0X)|2= G 2 F 2 ¯ `(k0)γµ(1−γ 5₎ 2 `(k)¯`(k)γ ν(1−γ 5₎ 2 `(k 0₎ (2.18) ×hhX|Jµ(q) |Pi hP|Jν†(q) |Xi i ,

where we identify the product of the expectation values of the electroweak currents as the hadronic tensor. It is then convenient to identify the term resulting from the leptonic vertex as the leptonic tensor Lµν_{, so that we can separate the leptonic and hadronic contributions}

to the cross section. The cross section is unpolarized, since the experiment studied in this project does not account for spin. Therefore the measured cross section, which is the cross section we are computing here, is an average over the incoming and outgoing lepton spins. Since we sum over all spin contributions we can write the leptonic tensor as

Lµν ₌ 1 4Tr h (k/0+m`)γµ(1−γ5)/k γν(1−γ5) i . (2.19)

We can then use the Dirac algebra and general properties of the trace to rewrite the leptonic tensor. This is done in the following way:

4Lµν ₌_Trh₍_k/0₊_m `)γµ(1−γ5)/k γν(1−γ5) i (2.20) =Tr /k0γµ/k γν−Tr /k0γµ/k γνγ5−Tr /k0γµγ5/k γν+Tr /k0γµγ5/k γνγ5 (2.21) +Tr[m`γµ/k γν] −Trm`γµ/k γνγ5 −Trm`γµγ5/k γν +Trm`γµγ5/k γνγ5 . (2.22)

Each term containing m`is zero, because these terms consist of an odd number of γ-matrices.

By inspection we see that the first and fourth term give equal contributions of the form Tr[γαk 0 αγ µ γβkβγ ν_{] =}_Tr_[ γαγµγβγν]k 0 αkβ (2.23) =4(gαµ_gβν₋_gαβ_gµν₊_gαν_gβµ₎_k0 αkβ (2.24) =4(k0µ_kν_{− (}_k0_·_k₎_gµν₊_k0ν_kµ₎_. _(2.25)

The second and third terms can be summed together and read

Tr /k0γµ/k γνγ5= −4ik0_αkβεαµβν. (2.26)

By collecting all contributions we find the full leptonic tensor: Lµν ₌₂₍_k0µ_kν₊_k0ν_kµ_{− (}_k0_·_k₎_gµν₋_ik0

αkβeαβµν). (2.27)

The final input needed to determine the square of the scattering amplitude is the hadronic tensor Wµν. As mentioned above, there are no Feynman rules for this object, thus we write

the hadronic tensor in terms of structure functions, F(x, Q2), that encode the partonic struc-ture of the nucleon. The hadronic tensor in its most general Lorentz covariant way reads [14] Wµν= −gµν+ qµqν q2 F1(x, Q2) + ˆ PµPˆν p·qF2(x, Q 2_{) −}_ie µνγδ qγ_pδ 2p·qF3(x, Q 2₎ _(2.28)

(12)

with ˆPµ= pµ−

p·q

q2 qµ.1 This parametrization is constructed by writing the hadronic tensor in

a linear independent fashion using the complete set of independent relevant tensors. We can decompose the hadronic tensor into a symmetric and anti-symmetric part: Wµν ₌ _Wµν

S +

W_Aµν. The relevant rank-2 tensors are gµν_{, q}µ_qν_{, p}µ_pν_{, q}µ_pν _{and the Levi–Civita tensor ε} αβµν.

We can write the symmetric part as a linear combination: W_Sµν = Agµν₊_Bqµ_qν₊_Cpµ_pν₊ D

2(p

µ_qν₊_pν_qµ₎_, _(2.29)

where A, B, C and D are arbitrary constants. Furthermore, we can add the anti-symmetric term due to the interference between vector and axial parts of the current [15]

W_Aµν =ieµνγδ

qγ_pδ

2p·qE, (2.30)

with E an arbitrary constant. By imposing current conservation, qµWµν =qνWµν =0, in the

hadronic vertex we can find the expressions for the arbitrary constants and write Wµν as in

Eq. (2.28). The full derivation can be found in appendix C. We can now find the scattering amplitude by contracting the leptonic tensor with the hadronic tensor:

Lµν_W µν= h 2(k0µ_kν₊_k0ν_kµ_{− (}_k0_·_k₎_gµν₋_ik0 αkβe αβµν₎i_× _(2.31) −gµν+ qµqν q2 F1(x, Q2) + ˆ PµPˆν p·qF2(x, Q 2_{) −}_ie µνγδ qγ_pδ 2p·qF3(x, Q 2₎ (2.32) The differential cross section can be simplified by a change of variables to DIS kinematics. After some algebra, which is written in appendix D, we find:

d2σνN dxdy = 4πα2 xyQ2ηW xy2F1(x, Q2) + (1−y− x2y2M2 Q2 )F2(x, Q 2_{) +}_y₍₁₋y 2)xF3(x, Q 2₎ , (2.33) where ηW is the contribution coming from the W+boson [14]. We have a simple expression

for the cross section with the structure functions that encode the information of the hadronic structure. The next step is to find an expression for these structure functions.

2.3 Parton Model

The Parton Model was proposed by Richard Feynman in 1969 [16] to describe high-energy hadron collisions. In the Parton Model it is assumed that the nucleon consists of different types of spin-1/2 point-like partons. The partons are then described by their momentum distributions. This means that the distribution of parton i is given by fi(ξ), where ξ is the

fraction of the longitudinal momentum of the parent hadron carried by the struck parton. The probability distributions, or parton distribution functions, are such that fi(ξ)dξ

repre-sents the number of i partons having a fraction between ξ and ξ+dξ of the nucleon longi-tudinal momentum. This analysis is performed in the infinite momentum frame (IMF), or

1_Pˆ

µ and pµ are different four momenta; the hat denotes that the vector describes an unobservable partonic

(13)

Breit frame, in which the proton has close to infinite longitudinal momentum, very small transverse momentum and the virtual boson has vanishing energy q0. This means we can write their respective four momenta as

pµ_{= (}_p0_{, 0, 0, p}

z)with pz ' p0(since we are in the deep inelastic regime) (2.34)

qµ_{= (}_q0 ₌_{0, 0, 0,}₋_q

z)with Q2 = −q2=q2z. (2.35)

In the IMF the interaction time of parton-parton interactions is time dilated, which means that the rate at which a parton interacts with its neighbouring parton during the time that the virtual boson is interacting with a parton is slowed down significantly. The parton is effectively a free particle at the instance of interacting with the boson. We can now assume the parton to be a free point-like particle, that it is carrying a fraction of the momentum of the parent hadron. The interaction of the boson with the parton is therefore elastic, so we can write:

ˆp2 = (ˆp+q)2 → Q2 =2 ˆp·q. (2.36) As the parton four-momentum is a fraction of the hadron four-momentum we can write it as ˆp =ξ p, which results in

2 ˆp·q=2ξ p·q= Q2 → ξ = Q

2

2p·q = x. (2.37) Hence the momentum fraction is identified as the Bjorken x as was mentioned earlier. Ac-cording to Feynman, in the Parton Model one can assume that the cross section is the inco-herent sum of partonic cross sections convoluted with the PDFs, so that we can write [13]:

dσ=

∑

i

d ˆσ fi(x)dx. (2.38)

where we replaced ξ with the Bjorken x variable. The partonic cross sections describe point-like particle interactions and can thus be computed using conventional methods. Eq. (2.38) shows that the hadronic cross section, which describes hadronic physics, can be written in terms of PDFs. Therefore we can write neutrino–nucleon scattering cross sections as [17]

d2_σνN

dxdy = G2_Fs

2π [x f(x) +x ¯f(x)(1−y)

2_]_, _(2.39)

This equation shows that we can express the hadronic cross sections not only in terms of structure functions, as in Eq. (2.33), but also in terms of the PDFs. We will see in section (2.4) how the structure functions are related to the PDFs.

2.3.1 Asymptotic freedom

The Parton Model uses the IMF to describe the quasi-free behaviour of the partons. This behaviour can however be more rigorously explained when introducing asymptotic free-dom, which is a characteristic property of Quantum chromodynamics (QCD). Asymptotic freedom ensures that the coupling αs= g2/4π, where g is the Yang-Mills coupling, between

(14)

partons becomes asymptotically weaker as the energy Q at which the hadron is probed in-creases. The β-function for a SU(3) symmetric system with nf massless fermions is written

as [17] β(g) = −g 3 (4π)2β0, with β0 =11− 2 3nf. (2.40) Then the equation at leading order in perturbation theory that shows the ‘running’ of the coupling reads

αs(Q2) =

αs(µ2)

1+β0αs(µ2)log(Q2/µ2), (2.41)

where µ is the renormalization scale. From this we can simplify the equation to [13]

αs(Q2) =

1

β0log Q2/Λ2_QCD

, (2.42)

whereΛQCD '300 MeV is the typical QCD scale. By inspection it is clear that αsis a function

of the momentum transfer Q and for Q2 6 ΛQCD the coupling αs > 1 and the partons are

tightly bound in the hadron. In the opposite scenario, if Q2 _Λ

QCD then αs 1 and the

partons are effectively free particles. The latter is referred to as asymptotic freedom and this is the limit where perturbative QCD calculations are valid. The first case is why the partons are unobserved in isolation and results in the color confinement property of QCD. [17]

2.4 QCD Improved Parton Model

The Parton Model described in 2.3 implies that the PDFs are effectively independent of Q2

in the DIS limit. Increasing the energy in which a hadron is probed implies the improvement of the spatial resolution, which means that at some energy the point-like substructure of the hadron is revealed. This means that increasing the energy does not further reveal any point-like structures, hence the dependency on Q2vanishes. This property is the so-called Bjorken scaling and is in fact only valid for leading order PDFs. It follows from DIS experimental data that this scaling is violated at small values of x, as is illustrated in figure 2.2. This means that the naive Parton Model is in fact not complete, as the transverse momentum is not restricted to be small. This is due to the parton being able to emit a gluon and therefore acquiring transverse momentum. The Parton Model must then be improved by QCD, by including the dynamics of the partons corresponding to the exchange of the gluons.

As we saw from the derivation of the cross section, the structure functions Fi(x, Q2)are

related to the PDFs inside the proton. Since the coupling is weak due to asymptotic free-dom, the structure functions can be separated by the use of the factorization theorem [18]. This means that the structure functions are written as a convolution of the PDFs with the coefficient functions. The coefficient functions C(x, αs(Q2))correspond to the perturbative

(15)

CHAPTER 2. DEEP INELASTIC SCATTERING reads [19] Fi(x, Q2) = nf

∑

i Ci(x, Q2) ⊗q(x, Q2) (2.43) = nf

∑

i x 1 Z x dy y Ci( x y, Q 2₎_f i(y, Q2),

where the sum runs over the different parton flavours and f is the density of a given parton. The explicit form of the structure functions is a result of the Mellin convolution, which is defined as [20] [f⊗g](x) = 1 Z x dz z f x z g(x) = 1 Z x dz z f(x)g x z. (2.44)

Even though the PDFs are non-pertubative objects, their Q evolution is determined by per-turbative QCD. This evolution is described by the DGLAP equations [21], which are given by [17] d d log Qfg(x, Q 2_{) =} αs(Q2) π 1 Z x dz z Pg←q(z) nf

∑

i fi( x z, Q 2_{) +} _¯f i( x z, Q 2₎ +Pg←g(z)fg( x z, Q 2₎ , (2.45) d d log Qfi(x, Q 2_{) =} αs(Q2) π 1 Z x dz z Pq←q(z)fi( x z, Q 2_{) +}_P q←g(z)fg( x z, Q 2₎ , (2.46) d d log Q ¯fi(x, Q 2_{) =} αs(Q2) π 1 Z x dz z Pq←q(z)¯fi( x z, Q 2_{) +}_P q←g(z)fg( x z, Q 2₎ , (2.47)

where f and ¯f are the parton and anti-parton distributions that can be treated as massless at scale Q and the sum runs over all quark flavours. The splitting function Pa←bis interpreted

as the probability density of a parton b emitting a parton a carrying momentum fraction xz. The explicit expressions for the splitting functions are

Pq←q(z) = 4 3 1+z2 (1−z)+ +3 2δ(1−z) , (2.48) Pg←q(z) = 4 3 1+ (1−z)2 z , (2.49) Pq←g(z) = 1 2[z 2_{+ (}₁₋_z₎2_]_, _(2.50) Pg←g(z) =6 1 −z z + z (1−z)+ + 11 12− nf 18 δ(1−z) . (2.51) The 1/(1−z)+term is defined in such a way that it agrees with the function 1/(1−z)for

(16)

the singular denominator 1/(1−z)that originally appears in the splitting functions in Eq. (2.48) and Eq. (2.51) [17]: 1 Z 0 dz f(z) (1−z)+ = 1 Z 0 dz f(z) − f(1) (1−z) . (2.52)

At this point we have written the DGLAP equations in the physical basis. This makes them complicated to solve, because all quarks couple to the gluon. It is convenient to decouple the non-singlet combinations of parton distributions from the gluon by writing the DGLAP equations in a different basis. We define the difference of parton distributions f_i± = fi± ¯fi,

where i is any type of quark. Then we can define the non-singlet distribution [20]:

f_NS,ij± (x, Q2) = (f_i±− f_j±)(x, Q2), (2.53)

where i, j denote any type of quark. The DGLAP equation for the non-singlet distribution reads d d log Q2fNS(x, Q 2_{) =} αs(Q2) 2π 1 Z x dz z fNS( x z, Q 2₎_P NS(z, Q2), (2.54)

where PNS is the non-singlet splitting function. The remaining distribution is called the

singlet distribution and is defined as the sum of all quark and antiquark flavours: Σ(x, Q2) =

nf

∑

i

f_i+(x, Q2) (2.55)

The singlet DGLAP equation is coupled to the gluon distribution g and is written as d d log Q2 Σ(x, Q2) g(x, Q2) ! = αs(Q 2₎ 2π 1 Z x dz z Pq←q Pq←q Pg←q Pg←g ! Σ(x_z, Q2) g(x_z, Q2) ! . (2.56)

(17)

Figure 2.2: The dependence on Q2of F2as measured by separate collaborations. The various

curves show the variations of F2in different regions of x [22].

2.5 Structure Functions

The observable measured by the NOMAD collaborations is the ratio of integrated cross sec-tions and will be described in detail in chapter 5. These cross secsec-tions are linear funcsec-tions

(18)

of structure functions and in this section we will define the relevant ones in this study. The structure functions can be split into the light, charm and heavy parts of the total proton struc-ture functions that read:

F₂ν,p= F_2,lν,p+F_2,cν,p+F_2,bν,p+F_2,tν,p (2.57) The light structure functions contain only parts that are proportional to the Cabibbo– Kobayashi–Maskawa (CKM) matrix elements corresponding to the light quark flavours; u, d and s quarks. The heavy part, including charm, is proportional to the CKM matrix elements Vh f or Vf h in such a way that mh > mf. The CKM matrix is a 3×3 unitary matrix which

contains the information on the strength of the quark-flavour changing weak interaction and relates the three quark generations.

The explicit expressions of the light and charm parts are [23] F_2,lν,p=2xC2,q⊗ |V_ud|2_d_{+ (|}_V ud| + |Vus|2)u¯+ |Vus|2s (2.58) +2 |Vud|2+ |Vus|2)C2,g⊗g , F_2,cν,p=2xC2,q⊗ |V_cd|2₍_d₊ _¯c_{) + |}_V cs|2(s+ ¯c) (2.59) +2 |Vcd|2+ |Vcs|2)C2,g⊗g ,

where C2,q are the perturbative coefficient functions. The remaining expressions are found

in [23]. Using the Callan-Gross relation it is straightforward to find the F1expressions. The

parity violating structure function F3is written in the same fashion [23]:

xF_3,lν,p =2xC3,q⊗ |V_ud|2_d_{− (|}_V ud| + |Vus|2)u¯+ |Vus|2s (2.60) +2 |V_ud|2_{+ |}_V us|2)C3,g⊗g xF_3,cν,p =2xC3,q⊗ |V_cd|2₍_d₋_¯c_{) + |}_V cs|2(s+¯c) (2.61) +2 |Vcd|2+ |Vcs|2)C2,g⊗g 2.5.1 Charm production

As explained in section 2.1 we neglect the masses of the partons in the initial and final states. However, in the case of charm production there is a charm quark in the final state, whose mass is not neglected. This means that if ˆp = ξ pis the parton four-momentum and there is

a charm quark in the final state, we find

m2_c = (ˆp2+q)2= /ˆp2−Q2+2 ˆp·q= −Q2+ξ Q

2

x (2.62)

This results in a different variable, the slow rescaling variable, that takes into account the en-ergy to excite the charm quark, on which the parton distributions above will be probed:

ξ =x 1+ m 2 c Q2 (2.63)

(19)

3. The NNPDF approach to parton

distribution functions

The PDFs are non-perturbative objects and must be extracted from experimental data. This means that the PDFs need to be parameterized to determine the x-dependence at an initial Q2₀, since the evolution in Q can then be determined by the DGLAP equations. The general form of the parametrization reads [24]

fi(x, Q20) = Nxαi(1−x)βiIi(x, a) (3.1)

where N is a normalization factor that ensures the sum rules on the PDFs. The power-like factors xαi _and₍₁₋_x₎βi _{describe the low-x and high-x behaviour of the PDFs and are}

theoretically motivated [25]. I (x, a) is a function that interpolates between the low-x and high-x regions. In general the interpolation function is chosen as some polynomial function dependent on a set of parameters{a}.

An alternative approach to determine the interpolation function is performed by the NNPDF collaboration that does not assume a polynomial form of the function but considers

I (x, a)to be the output of a neural network. This method reduces the bias from the choice of parametrization [25]. In this section we discuss two commonly used methods for error propagation: the Hessian method and the Monte Carlo method. In the case of the polyno-mial parametrization it is possible to use either the Hessian or the Monte Carlo approach to determine the uncertainties. For a neural network it is required to use the Monte Carlo method, because the Hessian approach is not affordable. This will become clear in section 3.2. Besides describing these two methods of error propagation, we give a basic summary of the neural network used by NNPDF and discuss the NNPDF3.1 PDF set that is the baseline of this project.

3.1 Monte Carlo method

The Monte Carlo method used by NNPDF to determine the uncertainties consists of the gen-eration of artificial datasets. This means that a set{F_i}of Nrepprobability distributions, i.e.

replicas, is generated using a Monte Carlo sampling method from the original experimental data. The replica generation is based on the experimental averages and information of the experimental covariance matrix. For a given experimental measurement the k-th replica is

(20)

CHAPTER 3. THE NNPDF APPROACH TO PARTON DISTRIBUTION FUNCTIONS

generated as [26]

F_i(art)(k)= S(_i,Nk)F_iexp

"

1+r_i(k)σ_istat+

Nsys

∑

j=1

r(_i,jk)σ_i,jsys

#

(3.2) with i = {1, . . . , Ndat}and j= {1, . . . , Nrep}, where Nrepand Ndatare the number of replicas

and number of values of the observable,F_i(exp), respectively. The normalization prefactor is defined as S(_i,Nk) ≡ N

∏

n=1 (1+r_i,n(k)σ_i,nN) (3.3)

The numbers r(_ik), r(_i,jk)and r(_i,nk)are randomly generated numbers from a univariate Gaussian distribution with the same correlation as the uncertainties [24]. At this point one has con-structed a Monte Carlo representation of the data, for which the mean value and variance of

F_i can simply be calculated by averaging over the replicas [27]:

E[F_i] = hF_ii. Var[F_i] = hF_i2i − hF_ii2 _(3.4)

The question is how many replicas are necessary to construct a faithful representation of the underlying probability density. It was found that with Nrep ∼ O(101)it is possible to

compute central values, Nrep ∼ O(102)allows to find uncertainties and Nrep ∼ O(103)is

required to properly reproduce data correlations [28] [29].

3.2 Hessian method

The Hessian method is a more conventional method for error propagation and both the ABM12 [30] and CT18 [31] PDF sets, which are fits from other collaborations, use versions of this method. A prediction of the observable using these sets is plotted in section 7.1.1 to make a comparison with the NNPDF predictions.

The Hessian method was first developed in [32] and is explained in detail in [33]. For our purposes we follow the discussion as in [4] and [29]. Given the χ2estimator, the best-fitting values of the parameters assume the χ2_{is quadratic about the global minimum. This means}

we can expand the estimator about the χ2_min: ∆χ2_≡ χ2−χ2_min = N

∑

i,j=1 H_ij(a_i−a0_i)(a_j−a0_j), (3.5)

where N is the number of fit parameters denoted by ai. The a0_i are the values of the fit

parameters that minimize the χ2. In Eq. (3.5) we have introduced the Hessian matrix Hij

whose components are the second derivatives of the χ2with respect to the fit parameters: Hij ≡ 1 2 ∂2χ2 ∂ai∂aj {a}={a0} (3.6) The uncertainties are then found, for an arbitrary functionF, by linear error propagation

σF =T N

∑

i,j ∂F ∂ai (H)−1∂F ∂aj !1/2 (3.7)

(21)

where T=p∆χ2_{is the tolerance factor for the required confidence interval that determines}

the range of allowed variations of the fit parameters and the 68% confidence interval. The inverse of the Hessian matrix is defined as the covariance matrix in the parameter space. To obtain it we need to invert the Hessian matrix, which is numerically not affordable when using a neural network since the dimensions of the matrix are significantly larger than the dimensions resulting from assuming a polynomial interpolation function.

The disadvantages of Eq. (3.7) are the derivatives of the observable with respect to the fit parameters, which are generally unknown. It is however possible to diagonalize the co-variance, or Hessian, matrix and represent the uncertainties in terms of a set of eigenvectors

{S±_i }[33]. After determining the eigenvectors, which is performed in more detail in [4], we find a simplified expression:

σF = 1 2 N

∑

i,j F (S+_i ) − F (S_i−) 2!1/2 , (3.8)

where S±_i indicate the i-th eigenvector with the positive or negative variation of the values minimizing the χ2.

3.3 Neural network

The neural network used by NNPDF to parameterize the PDFs is a multilayer feed-forward neural network, also known as a Perceptron [34]. The neural network provides unbiased determination of the interpolation function. It consists of a set of interconnected neurons ordered in input, hidden and output layers. Fig. 3.1 shows a schematic diagram of this ordering. The activation of the neurons is parameterized by a real number ξi, which in turn

is a function of the activation of neurons connected to the i-th neuron. Each pair of neurons is connected by a synapsis, which is assigned a sensitivity, or weight wij, with respect to the

output of the previous layer’s neuron. As the neural network is a feed-forward network, the layers are only influenced by the previous layers. The activation of a neuron in the(l+1)-th layer is then given by the recursive relation [35]

ξ_i(l+1)= g n_l

∑

j w(_ijl+1)ξ_j−θ_i(l) , with i= 1, . . . , nl, l=1, . . . , L, (3.9)

where θi is the activation threshold, or bias, of a given neuron i in layer l. This means there

is an input vector ξ1 _{in l} ₌ _{1 that is updated by every layer into the output vector ξ}L _with

l=L. We choose the activation function g to be a logistic function g(x) = 1 1+e−x, (3.10) where x= _∑nl j w (l+1) ij ξj−θ (l)

i . The nonlinearity of the sigmoid allows the neural network to

learn the nonlinear relations between the input and output and therefore produce nontrivial functions.

(22)

The idea of a neural network is for it to mimic the human brain, which means that it should be able to learn and adapt. There is therefore no need to know the underlying model that describes the data, which is exactly the reason why neural networks are useful. There are several types of training, or learning, algorithms that update the values of the weights and biases such that the network can reproduce an output. NNPDF uses a genetic algorithm to optimize the parameters in the neural network. A genetic algorithm works by repeatedly applying three operations [36]:

1. Selection, 2. Crossover, 3. Mutation.

The first step is to store the information of the parameter space into a chromosome. A point in this parameter space is then represented as

a= (a1, a2, . . . , anpar), (3.11)

where nparis the number of parameters. The bits ai correspond to the weights and biases of

the neural network. Once the chromosome is randomly initiated it is replicated ntot times.

Each chromosome is then assigned a f itness, which is the figure of merit that represents how close the given chromosome is to the best possible solution. The fitness is defined as the inverse of the χ2(a).

The goal is to find the best possible combination of parameters, which means that we repeatedly apply the aforementioned operations on the set of chromosomes. The first oper-ation is the mutoper-ation of a chromosome. A randomly selected bit is mutated by

ak →ak+η r− 1 2 , (3.12)

where the k-th bit is mutated, r is a uniform randomly generated number between 0 and 1, and η is the size of the mutation (mutation rate). The optimal size of the mutation rate can be dynamically determined and is dependent on the specific problem.

The second step is the crossover operation where random pairs of individuals are se-lected. For each pair, a random crossover bit s is chosen, and from this site the bits of the paired individuals are interchanged. The crossover step helps in obtaining a genetically dif-ferent child generation than the parent generation.

Finally, the chromosomes with best fitness are selected. This operation ensures that the best chromosomes propagate to the next generation. The procedure of subsequent opera-tions is an iterative process. Each iteration is called a generation and continues until the stopping criterion is satisfied. The stopping criterion is determined by means of a cross-validation method. The data is separated in two sets: the cross-validation set and the training set. The genetic algorithm is used to minimize the χ2 on the training set, while the χ2 on the validation set is monitored. The optimal set of parameters is achieved when the χ2 on the validation set stops improving [37].

(23)

Figure 3.1: Schematic overview of a multilayer feed-forward neural network.

3.4 NNPDF3.1

NNPDF3.1 is the set of PDFs used as the baseline in this project to compute all relevant quan-tities and is the set that will be improved after inclusion of the NOMAD data set. NNPDF3.1 is a global set of PDFs containing different kinds of data that are sensitive to the internal structure of hadrons, e.g. DIS, fixed-target Drell–Yan, W, Z and jet collider data. Addition-ally there is data gathered by ATLAS, CMS and LHCb from the LHC Run I [7] included in NNPDF3.1. Figure 3.2 shows the PDFs of the valence quarks, the strange and charm sea, and the gluon. The PDFs are shown at Q2 =10 GeV2on a logarithmic x-axis that shows the relatively large uncertainty of the strange-quark pdf.

Figure 3.2: The d−, u−, s+, c+quark and gluon g PDFs of NNPDF3.1 plotted at Q2=10 GeV2 in the range 10−4 ₆ _x

Bj 6 0.8. The plot clearly illustrates the relatively large uncertainty of

(24)

3.4.1 Data included in NNPDF3.1

There are several datasets included in the NNPDF3.1 fit, however only a few are sensitive to the strangeness in the same way as the NOMAD dataset is sensitive to the strangeness in the proton. The strangeness-sensitive datasets are tabulated in table 3.1. For the full table containing the datasets included in the NNPDF3.1 set, see [7]. The NuTeV dataset is the most similar to the NOMAD dataset, as both contain data from a DIS experiment. ATLAS is a collider experiment and as mentioned in chapter 1 the collider experiments show a significant difference in strangeness compared to DIS.

Experiment Obs. Ref. Nrep Kin1 range Kin2 range (GeV) Type

NuTeV σ_νcc [1][2] 45 (39/39) 0.026x60.33 2.06 Q610.8 DIS σcc_¯ν [1][2] 45 (37/37) 0.026x60.21 1.96x68.3

ATLAS W, Z 2010 [38] 30 (30/30) 06 |ηl| 63.2 Q= MW, MZ Collider

W, Z 2011 [9] 34 (34/34) 06 |η_l| 62.3 Q= MW, MZ

Table 3.1: Table showing the data sensitive to the strangeness in the proton included in the NNPDF3.1 fit. [7]

(25)

4. Bayesian reweighting and

unweighting

In general whenever new datasets are obtained from deep inelastic and hadronic collision experiments, it is possible to determine the PDFs and their uncertainties by fitting the new data. Including new data by fitting is in general a time consuming process. An alternative method to include new data to an ensemble of Nrep replicas of PDFs is by reweighting and

unweighting. By computing the χ2_{between the new data and the theoretical predictions of}

the new data, one can assign weights to each member of the ensemble of Nrep PDFs. The

weight is determined by the value of the χ2, which means that if the theoretical prediction agrees well with the new data the corresponding member of the ensemble is assigned a rel-atively large weight. Predictions that show poor agreement correspond to weights that go to zero. This results in a weighted ensemble that gives an updated distribution of the pre-diction of the new data, which is similar to fitting to new data. Reweighting is less time consuming than performing a new global fit, however it is not without consequences. The-oretical predictions that are assigned a small weight, which thereby agree poorly with the data, are effectively removed from the ensemble which then results in a loss of information. Following reweighting there is the unweighting procedure: the transformation of weights into the number of copies of a specific member of the original ensemble. This section ex-plains the reweighting and unweighting procedure and is largely based on Refs. [10] and [11] respectively.

4.1 Reweighting

First we start with a finite ensemble of PDFs ξ = {fk|k = 1, . . . , Nrep}, which corresponds

to the probability distribution of PDFs,P_prior(f), based on an arbitrary set of experimental data. It is possible, using the ensemble, to compute any statistical quantity or observable

O[ f ] depending on the PDFs. To update the prior when new experimental data become available, one can either perform a fit on all available data, or use reweighting and compute a weight, wk, for each fkin ξ depending on the agreement of the theoretical prediction with

the new data. All replicas in the original ensemble ξ start with wk = 1 implying all are

(26)

CHAPTER 4. BAYESIAN REWEIGHTING AND UNWEIGHTING

Nrep replicas, which then allows the calculation of the expectation value

hOi = Z O[f]P_prior(f)d f = 1 Nrep Nrep

∑

k=1 O[fk], (4.1)

where the integral is over the space of PDFs [10]. Considering a set of new data one wants to find the same expression as in Eq. (4.1) but with wk 6=1. For this, the weights need to be

calculated, which is done as follows. First it is assumed that the new dataset is statistically independent of the already included data. Then, assuming the dataset has n datapoints and systematic uncertainties are not correlated, there is an n×n experimental covariance matrix

σ_ij that describes the systematic uncertainties. By using statistical inference it is the goal to

update the prior and obtain an improved probability density function, the posterior. To do this, one needs to use the existing theory to compute the relative probabilities of the new data for each replica. Since the data are assumed to have Gaussian errors, these probabilities are directly proportional to the probability density of χ2given a replica f :

P (χ2|f)∝(χ2) 1 2(n−1)e−12χ2 (4.2) where χ2is evaluated by [29] χ2(f) = n

∑

i,j=1 (Di−Ti[f])σ_ij−1(Dj−Tj[f]) (4.3)

Here Di are the new datapoints and Ti[f]are the predicted values of the datapoints using f .

The covariance matrix has the general form

σij ≡σiσjδij+ n

∑

i,j=1

σiσj (4.4)

where δijis the Kronecker delta and the off-diagonal term is zero due to the absence of

corre-lated systematic uncertainties. We can apply the rules for multiplication of probabilities due to the statistical independence of the old and new data [10] and write the new probability distribution as

P_posterior(f) = N P (χ|f)Pprior(f) (4.5)

whereN is as of yet an undefined normalization.

Using the knowledge of the expression for the new probability distribution we compute the expectation value of the observable using the posterior:

hOi_new= Z O[f]P_posterior(f)d f = N Z O[f]P (χ|f)Pprior(f)d f = 1 Nrep Nrep

∑

k=1 N P (χ|fk)O[fk] (4.6)

(27)

where in the second line we replaceP_posteriorwith Eq. (4.5) and in the final line we used Eq. (4.1). From here it follows that we can compute the new expectation value as

hOinew= 1 Nrep Nrep

∑

k=1 wkO[fk], (4.7) such that w_k = N P (χ|f_k) = N0(χ2_k) 1 2(n−1)e−12χ2k _(4.8)

where χ2_k = χ2(fk)and corresponds to the χ2 between the prediction using fk and the new

datapoint. The normalization factorN0 _{is defined by considering the observable}_O[_f_]_{to be}

the unit operator which then imposes the condition Nrep =

Nrep

∑

k=1

wk (4.9)

This shows that the sum of the weights must equal the number of replicas. Using Eq. (4.8) and Eq. (4.9) we find the expression for the weights:

wk = (χ2_k)12(n−1)e−12χ2k 1 Nrep ∑ Nrep k=1(χ2k) 1 2(n−1)e−12χ2k (4.10)

Dividing this expression by Nrep gives us the probabilities of the replicas fk.

4.2 Loss of information

The maximally efficient representation of the underlying probability density P_prior is the original ensemble of replicas ξ where each replica has equal weight. This density can simply be improved by increasing the number of replicas Nrep. This is not the case after reweighting,

due to several replicas being effectively irrelevant in ensemble averages due to their small weights. There is therefore a loss of information, which means the accuracy of the new probability density is smaller than it would be if there would have been a completely new fit. The loss of information can be quantified by using the Shannon entropy which allows us to compute the number of effective replicas remaining after the reweighting procedure [10]:

Neff ≡exp " 1 Nrep Nrep

∑

k=1 wkln Nrep wk # (4.11) The number of effective replicas Neff is the number of replicas Nrep needed in a new fit to

obtain equal accuracy as the reweighted set. This means that if Neffis too low, reweighting

becomes unreliable. Neffdecreases if the new data contains a lot of new information on the

PDFs or if new data are inconsistent with the old data. If there is no new information and there is consistency with the old data not much is gained and Neff would be approximately

equal to Nrep. This means that the prior has to be large enough to result in a statistically valid

(28)

It is possible to determine whether the new data is inconsistent or if the new data contains a lot of information by computing the α-profile of the data. We define inconsistent data as data whose errors have been underestimated. These errors can be rescaled by a factor α. This rescaling is done by the transformation χ2

k → χ2k/α2, which results in a probability density

of α P (α)∝ 1 α Nrep

∑

k=1 wk(α) (4.12)

If this probability has a maximum around α = 1 the data are considered consistent and if it peaks far from one it is more probable that the data have underestimated uncertainties.

4.3 Unweighting procedure

The step following reweighting is to unweight the reweighted set so that it can be used without the inclusion of the weights for each different replica. We start with a set ξ of Nrep

reweighted replicas, that each carry a weight w_k as defined in Eq. (4.10). The aim of un-weighting is to obtain a new set of PDFs containing N_rep0 6= Nrep with all weights equal

to one. This new set will have the same probability distribution as the weighted set and will not provide new physical insights. It does however have practical purposes, since the unweighted set can then be used similarly as a set before reweighting. This means that reweighting and unweighting is in principle an endless cycle as long as there are enough replicas surviving.

The starting point for unweighting is again the set ξ of Nrep reweighted replicas. To

build the unweighted set we select replicas from ξ in such a way that replicas with higher weights are selected more often, while replicas with vanishing weights are not included in the unweighted set. Note that the choice of N_rep0 , the size of the unweighted set, is in principle arbitrary, the weighted set does not contain more information than N_effreplicas. This means that if N_rep0 > Neff results in a redundant set, as there is no new information included. We

therefore choose N_rep0 = Neff so there is no loss of information and no redundancy. The

procedure then continues with determining the probability of each replica pk =

wk

Nrep

, (4.13)

which follows from the normalization as seen in Eq. (4.9). The next step is to define proba-bility cumulants P_k ≡P_k−1+pk = k

∑

j=0 pj, (4.14)

where we take P0 = 0. It follows that 06 Pk 6 1 and Pk−1 6 Pk. We then determine the set

of Nrepnew weights as follows:

w0_k = N0 rep

∑

j=1 θ j N0 rep −Pk−1 ! θ Pk− j N0 rep ! (4.15)

(29)

where θ is the Heaviside step function that causes each new weight w0_kto be either a positive integer or zero. Furthermore, the new weights satisfy the normalization condition

N_rep0 =

Nrep

∑

k=1

w0_k. (4.16)

The unweighted set is then constructed by selecting w0_k copies of the k-th replica for all k ∈ {1, . . . , Nrep}. The probability of replica k in the unweighted set is now given by

p0_k = w 0 k N0 rep . (4.17)

We find the original probability of the weighted set by taking the limit N_rep0 → ∞ which

(30)

5. Charm production in NOMAD

The NOMAD collaboration has performed an experiment measuring the cross sections re-lated to neutrino–iron interactions resulting in charm dimuon production. After correcting for background data, they observed 15,344 charm dimuon events, which is the largest data sample currently available. The observable measured by the NOMAD collaboration, which is the central object in this study, will be described in this section. This discussion and the definition of quantities is mostly based on the NOMAD paper regarding this experiment. [39]

5.1 Charm dimuon production

The process of charm dimuon production is initiated by a νµ induced CC production of a

charm quark, which in turn semi-leptonically decays into a final state secondary muon µ+ with an opposite charge to the final state muon µ− of the leptonic vertex. Fig. 5.1 shows a Feynman diagram of this interaction. We see that the incoming quarks can be either strange or down quarks. Although there might be a down quark in the initial state, this contribution is suppressed by the CKM matrix elements. This means that the most significant contribu-tion to the cross seccontribu-tion σµµof this process is coming from the strange sea of the nucleon.

Figure 5.1: Feynman diagram corresponding to a νµinduced CC charm dimuon production

(31)

CHAPTER 5. CHARM PRODUCTION IN NOMAD

5.2 The observable: R

µµ

The NOMAD observable is defined as the ratio of the charm dimuon cross section, σµµ, with

the inclusive DIS neutrino–nucleon cross section, σCC:

Rµµ(x) ≡ σµµ σCC , (5.1) where x = Eν, xBj, √

ˆs and the partonic centre-of-mass squared is defined as ˆs = Q2₍₁₋

xBj)/xBj. In this thesis we focus solely on the dependency of the observable on the neutrino

energy Eν. The analysis on the other dependencies is beyond the scope of this thesis and left

for further study. As mentioned above, the cross section σµµis highly sensitive to the strange

content. This is in contrast to the inclusive cross section σCC which is relatively insensitive

to the strange content of the nucleon due to the process W−d→u not being suppressed in a similar way as in the charm dimuon production [40]. This means that the ratio defined in Eq. (5.1) is therefore a good measure for the strange content of the nucleon. Another advantage coming from the use of this ratio is the large cancellation of the systematic uncertainties affecting both the numerator and the denominator equally.

5.2.1 The denominator: σCC

The denominator is the inclusive neutrino induced CC cross section σCC. We give the

expres-sion of the differential cross section, which reads [12] d2σ_CCνN dxdQ2 = G2_FM4_W 4πx(Q2₊_M2 W)2 Y+− 2M2_Nx2y2 Q2 FνN 2,CC(x, Q2) +Y−xF3,CCνN (5.2) −y2FνN L,CC ,

where MW is the W boson’s mass and Y±=1± (1−y)2. The longitudinal structure function

FLis defined by the Callan–Gross relation in the improved Parton Model by FL =F2−2xF1.

To find the expression in Eq. (5.2) from Eq. (2.33) we use this extended Callan–Gross relation and a change of variables.

The structure functions used to compute the cross section are the total structure functions, which are defined in section 2.5. Note that the cross sections in the observable, and there-fore the structure functions, describe neutrino–nucleon interactions. The reason why the structure functions are not explicitly describing neutrino–iron interactions is because iron is approximately isoscalar. Due to SU(2) isospin symmetry relating neutrons and protons by the transformation u → d, all isoscalar targets give equal contributions. We take a nucleus structure function F2= 1 A ZF₂p+ (A−Z)F₂n , (5.3)

where F₂p(n)are the proton(neutron)’s structure functions, A is the atomic mass number and Z is the atomic number. When considering an isoscalar target we have Z = A−Z, which

(32)

CHAPTER 5. CHARM PRODUCTION IN NOMAD means A=2Z. Hence, F2= 1 2(F p 2 +F2n) (5.4)

for all isoscalar targets.

5.2.2 The numerator: σµµ

The differential cross section for charm production in CC neutrino DIS is given by [39] d2σ_cνN dxdy = G_F2MNEν π(1+Q2/M_W2 )2 (1−y− Mxy 2Eν Fν (2,c)(x, Q2) + y2 2 F ν T,c(x, Q2) (5.5) +y(1−y 2)xF ν 3,c(x, Q2) , where we use the charm structure functions as defined in section 2.5. In order to write this expression in a similar fashion as Eq. (5.2), we use the relation FT = 2xF1 [41], the DIS

kinematics and perform a change of variables to find d2σ_cνN dxdQ2 = G2_FM_W4 4πx(Q2₊_M2 W)2 Y+− 2M2_Nx2y2 Q2 FνN 2,c(x, Q2) +Y−xF3,cνN (5.6) −y2FνN L,c

In order to define the charm dimuon cross section which appears as the numerator, we need to introduce the semi-leptonic branching ratio Br(h →µX)and the fragmentation function

Dh_c of the charm quark into a given charmed hadron. The branching ratio corresponds to the fraction of decay modes from hadron h to final state µX. The branching ratio can be written implicitly as Bµ =∑h fhBr(h →µX), where fhis the production fraction of h, and is

a function of the neutrino energy. The explicit form is parametrized as Bµ(Eν) =

a 1+b/Eν

(5.7) where the parameters are fitted by NOMAD [39] to E531 data [42]. The values of the param-eters are a=0.097 ± 0.003 and b =6.7 ± 1.8 GeV.

The charm fragmentation function Dh_c(z)describes the probability that the charm quark fragments into hadron h with h carrying fraction z of the parent parton momentum. The ex-act form of the fragmentation function can be expressed by the Collins–Spiller parametriza-tion [43]: Dc(z) = 1 −z z +εc 2−z 1−z (1+z)2 1−1 z − εc 1−z −2 (5.8) with εcbeing a free parameter. The value of this fragmentation function is however included

in the data published by NOMAD. We can therefore simplify and write the dimuon differ-ential cross section as

d2_σ µµ dxdQ2 =Bµ(Eν) d2_σνN c dxdQ2. (5.9)

(33)

6. Calculation of neutrino DIS cross

sections

In order to properly determine the impact of the NOMAD data on the PDFs one needs to re-produce the experimental dataset given by NOMAD. This dataset consists of the observable Rµµas a function of Eνas described in chapter 5. The theoretical calculation of the observable

is implemented in a standalone PYTHON program2using the APFEL library [44] to evaluate the structure functions. This section explains the benchmark done to verify the validity of the code, the expressions implemented in the code and the kinematic limits used for integration. Finally a brief summary of the APFEL library and the FONLL method, that includes heavy quark mass contributions to the structure functions, is given.

6.1 Benchmarking the code

The PYTHON program’s structure is based on the code that was used to compute the neutrino induced CC cross sections given in table A.1 of [12]. The upper window in figure 6.1 shows the distribution of the numbers as shown in table A.1 of [12] and the results of the code. The lower window shows the discrepancy between the numbers and the calculation. It is clear that the calculation done with the PYTHON code shows excellent agreement in the lower neu-trino energy Eν range, but deviates in the high energy range up to 5–10%. The benchmark

is therefore successful as the kinematic range of the NOMAD dataset lies in the low energy region. The discrepancy at ultra-high energies will then have negligible effect on the results. The results of the benchmark are obtained by using the FONLL scheme in NNLO with a maximum of nf =5 active flavours, and using the NNPDF3.1sx+LHCb NNLO+NLLx PDF

set. The original numbers were produced with nf = 6 active flavours, however the impact

of the t-quark is small and negligible for the kinematic range of the NOMAD dataset.

(34)

CHAPTER 6. CALCULATION OF NEUTRINO DIS CROSS SECTIONS

Figure 6.1: Upper: Plot showing the distribution of the numbers shown in table A.1 of [12] and the computation done with the PYTHON program. Lower: Plot illustrates the discrep-ancy between the theoretical calculation and the numbers computed in [12]. There is a clear agreement in the low energy range and the theoretical calculation shows a discrepancy of 5–10% at high energies.

6.2 Calculation of the cross section

6.2.1 The integrands

As explained in chapter 5 the NOMAD observable Rµµ is the ratio of the semi-inclusive

dimuon production neutrino–nucleon cross section with the inclusive CC cross section. Both cross sections are computed separately in two similar PYTHON programs by integrating their respective expressions. The integrand reads

d2_σ CC,c(x, Q2, Eν) dxdQ2 = G2_FM_W4 4πx(Q2₊_M2 W)2 Y+− 2M2_Nx2y2 Q2 FνN 2,(CC,c)(x, Q2) +Y−xF3,νN(CC,c) (6.1) −y2FνN L,(CC,c)

(35)

CHAPTER 6. CALCULATION OF NEUTRINO DIS CROSS SECTIONS

which corresponds to the equation derived in Eq. (2.33). This expression holds both for inclusive scattering and charm production scattering. The difference is the form of the struc-ture functions as is explained in section 2.5. The charm production cross section contains only the Fc structure functions and the inclusive cross section is computed using the total

structure functions. However, to obtain the dimuon cross section, which is the numerator of Rµµ, one needs to include the branching ratio Bµ in the charm production calculation as

defined in Eq. (5.7):

d2σµµ

dxdy =Bµ d2σc

dxdy. (6.2)

The integral is computed as a function of Eνwith the values for the energy given by

NO-MAD in table 6 of [39], which is tabulated in table B.1. The structure functions are evaluated using APFEL, which will be explained in section 6.2.3. The numerical values of the heavy quark masses implemented in the code depend on the PDF set that is used in the computa-tion. However, GF= 1.663787·10−5 GeV−2[14] is constant in all calculations. To construct

the theoretical predictions we integrate over the allowed kinematic range in x and Q2. The integral then reads

σ(Eν) = Z Q2_max Q2 min dQ2 Z x_max xmin dx d 2 σ dxdQ2(x, Q 2_{, E} ν) (6.3) This is evaluated numerically using adaptive quadrature as it is implemented in the SciPy library with the relative precision set to 10−8. Note that (6.3) is a one-dimensional nested integral. We use this expression because it results in an improved computational efficiency compared to evaluating the double integral. To further improve the efficiency one needs to perform a change of variables

x→log x, Q2 →log Q2 (6.4) This results in the integral

σ(Eν) = Z log Q2_max log Q2 min d log Q2 Q2 Z log x_max log xmin d log x x d 2 σ dxdQ2(x, Q 2_{, E} ν) (6.5) 6.2.2 Kinematic limits

It is important to notice that the integrated variables are arbitrarily chosen and using the DIS kinematics it is possible to compute the integral in any other combination of the DIS variables. The change of variables reads

d2_σCC νN(x, Q 2_{, E} ν) dxdQ2 = 1 xs d2_σCC νN(x, y, Eν) dxdy = y x d2_σCC νN(Q 2_{, y, E} ν) dQ2_dy (6.6)

where the relations of the DIS kinematics are used as explained in section 2.1. As seen in Eq. (6.1) the cross section depends on x, Q2and Eνand from the definition of the inelasticity, y,

we know 06y61 and the relation y= Q_xs2 3_{. The constraints on y impose relations between}

the variables and determine the integration limits in the following way: 3_s₌_2M

Impact of the NOMAD data on the strangeness in the proton

T

HESIS

MS

C

P

ROJECT

T

HEORETICAL

P

HYSICS

Impact of the NOMAD data on the

strangeness in the proton

BY

FERRAN

FAURA

I

NATIONAL

I

S

P

Contents

1. Introduction

2. Deep Inelastic Scattering

2.1

DIS Kinematics

2.2

Differential cross section

2.3

Parton Model

∑

2.4

QCD Improved Parton Model

∑

∑

∑

∑

2.5

Structure Functions

3. The NNPDF approach to parton

distribution functions

3.1

Monte Carlo method

∑

∏

3.2

Hessian method

∑

∑

∑

3.3

Neural network

∑

3.4

NNPDF3.1

4. Bayesian reweighting and

unweighting

4.1

Reweighting

∑

∑

∑

∑

∑

∑

4.2

Loss of information

∑

∑

4.3

Unweighting procedure

∑

∑

∑

5. Charm production in NOMAD

5.1

Charm dimuon production

5.2

The observable: R

6. Calculation of neutrino DIS cross