Measurement of the W boson mass

The $W$ boson mass is measured using proton-proton collision data at $\sqrt{s}=13$ TeV corresponding to an integrated luminosity of 1.7 fb$^{-1}$ recorded during 2016 by the LHCb experiment. With a simultaneous fit of the muon $q/p_T$ distribution of a sample of $W \to \mu\nu$ decays and the $\phi^*$ distribution of a sample of $Z\to\mu\mu$ decays the $W$ boson mass is determined to be \begin{equation*} m_{W} = 80354 \pm 23_{\rm stat} \pm 10_{\rm exp} \pm 17_{\rm theory} \pm 9_{\rm PDF}~\mathrm{MeV}, \end{equation*} where uncertainties correspond to contributions from statistical, experimental systematic, theoretical and parton distribution function sources. This is an average of results based on three recent global parton distribution function sets. The measurement agrees well with the prediction of the global electroweak fit and with previous measurements.


Introduction
The W boson mass (m W ) is directly related to electroweak (EW) symmetry breaking in the Standard Model (SM) [1][2][3].At tree level, m W = gv/2 where g is the weak-isospin coupling and v is the vacuum expectation value of the Higgs field.Going beyond tree level the boson masses and couplings receive loop corrections.The value of m W is related to the precisely measured fine-structure constant (α), the mass of the Z boson (m Z ) and the Fermi constant (G F ), as [4,5] where ∆ encapsulates the loop-level corrections.A global fit of EW observables, excluding direct measurements of m W , yields a prediction of m W = 80354 ± 7 MeV [6]. 1 This can be compared with direct measurements to test for possible beyond SM contributions to ∆ in Eq. 1.The 2020 PDG average of direct measurements is m W = 80379 ± 12 MeV [7].
The sensitivity of the global EW fit to physics beyond the SM is primarily limited by the precision of the direct measurements of m W [6]. Furthermore, the uncertainty in the prediction is expected to reduce as the top-quark mass, which is the leading source of the parametric uncertainty, is determined more precisely in the future.The value of m W was measured to a precision of 33 MeV at the Large Electron-Positron (LEP) collider [8] at CERN and to a precision of 16 MeV in an average [9] of measurements by the CDF [10] and D0 [11] experiments at the Fermilab Tevatron collider.The first measurement at the LHC was performed by the ATLAS collaboration and has an uncertainty of 19 MeV [12].The hadron collider measurements are based on three observables in leptonic W boson decays, namely the transverse mass, missing transverse momentum and charged lepton transverse momentum (p T ).At hadron colliders, the lepton p T is measured with good resolution but it is strongly influenced by the W boson transverse momentum distribution, the modelling of which is a potential source of a limiting systematic uncertainty.However, the resolution of the transverse mass is degraded by the pile-up of proton-proton interactions in the same bunch crossing.Therefore, the lepton p T was the most sensitive observable in the recent measurement performed by the ATLAS collaboration.Despite being based on a small subset of the data recorded to date, the ATLAS measurement of m W is already limited by uncertainties in modelling W boson production, in particular the parton distribution functions (PDFs) of the proton.
The potential for a measurement based on the muon p T with the LHCb experiment is studied in Ref. [13] It was estimated that LHCb data collected in LHC Run 2, at a proton-proton (pp) centre of mass energy √ s = 13 TeV, would allow a measurement with a statistical precision of around 10 MeV.Owing to the complementary pseudorapidity (η) coverage of the LHCb experiment with respect to the ATLAS and CMS experiments, it was demonstrated in Ref. [13] that the PDF uncertainty could partially cancel in an average of m W measurements by the LHC experiments.
In this paper a first measurement of m W is presented using W → µν decays, including both W boson and muon charges, collected at the LHCb experiment. 2 This measurement considers the muon q/p T distribution, where q is the muon charge.Figure 1 (left) illustrates 1 Throughout this paper natural units with c = 1 are used. 2The inclusion of charge-conjugate processes is implied throughout unless otherwise specified.
how the shape of the muon q/p T distribution in simulated W boson decays is influenced by variations in m W of ±300 MeV, which corresponds to roughly ten times the target precision of the present analysis.The q/p T variable allows all muons with p T > 24 GeV to be visualised; those with 28 < p T < 52 GeV are used to determine m W , while consistent control of the fit can be demonstrated in the region p T > 52 GeV.
The p T of a muon produced by the decay of a W boson has a strong dependence on the W boson transverse momentum (p W T ).Direct measurements of the p W T distribution have been reported by the ATLAS [14] and CMS [15] collaborations but the intervals are necessarily coarse due to the limited p W T resolution.Measurements of the transverse momentum distribution for Z boson production (p Z T ) are therefore used to validate the predictions for the p W T distribution. 3The angular variable φ * [16], defined in Eq. 4, is used in this analysis as a proxy for p Z T since its distribution can be measured more precisely than that of p Z T .Parton-shower programs such as Pythia [17] can be tuned (e.g.Ref. [18]) to describe the p Z T and φ * data at the per cent level but it is challenging to reliably translate such tunes to W boson production.However, a W -boson-specific tuning of a parton-shower model can be performed simultaneously with a determination of m W [19].
If electroweak corrections are neglected then the production and leptonic decay of the W boson factorise such that the differential cross-section can be written as dσ dp W  T dydM d cos ϑdϕ sin 2 ϑ cos 2ϕ + A 3 sin ϑ cos ϕ + A 4 cos ϑ + A 5 sin 2 ϑ sin 2ϕ + A 6 sin 2ϑ sin ϕ + A 7 sin ϑ sin ϕ , where ϑ and ϕ are the lepton decay angles defined in a suitable frame (the Collins-Soper frame [20] is used in this analysis), and p W T , y and M denote the transverse momentum, rapidity and mass of the final state lepton pair, respectively.An equivalent expression applies to Z → µµ production.The eight angular coefficients (A i ) are ratios of helicity cross-sections and depend on p W T , y and M ; σ unpol. is usually referred to as the unpolarised cross-section.The coefficients A 5 − A 7 are numerically small because they only arise at second order, or higher, in the strong coupling constant (α s ). 4 The coefficient A 3 is particularly influential on the muon p T distribution.Figure 1 (right) shows how the q/p T distribution in simulated W boson events, after the selection requirements described in Sect.3, changes when A 3 is scaled up and down by 3%.
In this paper, the simulated samples are weighted in the full five-dimensional phase space of vector boson production and decay, using different models for the unpolarised cross-section, angular coefficients, and QED final-state radiation.Several PDF sets are used in the analysis but none of the analysed data were included in the determination of these PDF sets.
This paper is organised as follows.Section 2 describes the data and simulated samples.Section 3 details the signal candidate selection requirements.Section 4 describes charge-dependent curvature corrections that are applied to the data and simulation.The determination of residual smearing corrections to the simulation with a simultaneous fit of Z → µµ and quarkonia decays is subsequently described.Section 5 details the measurement of muon selection efficiencies and subsequent weight-based corrections to the simulation.Section 6 describes the treatment of background arising from in-flight decays of light hadrons.Section 7 sets out the modelling of vector boson production and decay.Section 8 describes the simultaneous fit of the model to the muon p T distribution of W boson candidates and the φ * (defined in Eq. 4) distribution of Z boson candidates to determine m W . Section 9 explains how results based on three different PDF sets are averaged, and summarises the systematic uncertainties.Several cross-checks of the measurement are reported.The impact of analysis choices and systematic variations on m W is discussed throughout the paper.The conclusions of the analysis are presented in Sect.10.

Data sets and event selection
The LHCb detector [21,22] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks.The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the proton-proton (pp) interaction region [23], a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes [24] placed downstream of the magnet.The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV.The minimum distance of a track to a primary pp collision vertex (PV), the impact parameter, is measured with a resolution of (15 ⊕ 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV.Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors [25].Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter.Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers [26].The online event selection is performed by a trigger [27], which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction.
This analysis uses a data sample of pp collisions at √ s = 13 TeV recorded during 2016, corresponding to an integrated luminosity of about 1.7 fb −1 .Roughly half of the data were recorded in each of the dipole magnet polarity configurations, resulting in a large degree of cancellation of charge-dependent curvature biases and their associated uncertainties.These data correspond to an average number of proton-proton interactions per bunch-crossing event of O(1).
During Run 2 the LHCb detector was aligned and calibrated in real-time [28].The alignment of the tracking system is based on a χ 2 minimisation of the residuals of the clusters of tracker hits evaluated with a Kalman filter that takes into account multiple scattering and energy loss [29].The alignment algorithm also permits mass and vertex constraints [30].An optimised offline alignment, which includes Z → µµ events with a mass constraint that accounts for the natural width of the Z boson, is used to determine the track parameters.Since the real-time alignment is not optimised for the analysis of high-p T final states, this realignment improves the Z → µµ mass resolution by around 30%.
All signal and background processes are simulated using an LHCb specific tune [31] of Pythia version 8.186 [17].The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [32] as described in Ref. [33].Events are simulated with both polarity configurations and weights are assigned to events in each polarity such that the polarity distribution matches the recorded data.

Selection of W boson, Z boson and quarkonia signal candidates
Tracks are identified as muons if they are matched to hits in either three or all four of the most downstream muon stations depending on their momentum.They are then considered in this analysis if they are within the range 1.7 < η < 5.0, and have a momentum of less than 2 TeV.The tracks must have a good fit quality and a relative momentum uncertainty of less than 6%.Candidate W → µν events are selected online by requiring that one identified muon satisfies the requirements of all stages of the trigger.At the hardware stage a p T of at least 6 GeV is required.The isolation of a muon is defined as the scalar sum of the transverse momenta of all charged and neutral particles, as selected by a particle-flow algorithm, described in Ref. [34], within (∆η) 2 + (∆φ) 2 < 0.4 2 around the muon, where ∆η and ∆φ denote the separation in η and azimuthal angle around the beam direction (φ), respectively.Hadronic background contributions are suppressed by requiring the muon to have an isolation of less than 4 GeV.For the W → µν selection the η range is tightened to 2.2 < η < 4.4 so that the area of the isolation cone is fully instrumented.The muon must satisfy χ 2 IP < 9 where χ 2 IP is defined by the difference in the vertex fit χ 2 of the PV with and without including the muon.Background from Z boson events is suppressed by rejecting events that contain a second muon with p T > 25 GeV and an opposite charge to that of the primary muon candidate.Roughly 2.4 million W → µν candidates are selected in the range 28 < p T < 52 GeV.
Candidate Z → µµ events are reconstructed from combinations of two oppositely charged identified muons associated to the same PV with an invariant mass within ±14 GeV of the known Z boson mass [7].At least one muon must be matched to a single muon selection at all stages of the trigger.Both muons must have p T > 20 GeV, an isolation value below 10 GeV, and an impact parameter significance of less than ten standard deviations.Roughly 190 thousand Z → µµ candidates are selected.
Candidate J/ψ → µµ and Υ (1S) → µµ events, which are primarily used to calibrate the modelling of the momentum measurement, are required to have a pair of oppositely charged identified muons.Both muons must have a transverse momentum above 3 GeV and satisfy a tighter muon identification requirement.In order to specifically select J/ψ → µµ candidates originating from b-hadron decays the decay vertices must be displaced from the nearest PV with a significance of at least three standard deviations.These selections retain roughly 1.0 million Υ (1S) → µµ candidates and 220 thousand J/ψ → µµ candidates.

Momentum calibration and modelling
The momentum scale can be precisely determined from the mass measurements of various resonances, including those that decay to muon pairs.However, charge-dependent curvature biases that shift q/p are challenging to estimate because their effect largely cancels in the mass of the resonances.They are also particularly important for the high momentum muons from W and Z boson decays.In Ref. [35] it was proposed to determine corrections using the so-called pseudomass variable in Z → µµ events where p ± and p ± T are the momenta and transverse momenta of the µ ± , respectively.The opening angle between the two muons is denoted θ.Crucially, the value of M ± is independent of the magnitude of the momentum of the µ ∓ and is therefore directly sensitive to curvature biases affecting the µ ± candidate.The pseudomass is an approximation of the dimuon mass under the assumption that the dimuon system has zero momentum transverse to the bisector of the two lepton transverse momenta.The φ * observable is defined as [16] where ∆φ is the azimuthal opening angle between the two leptons and ∆η is the difference between the pseudorapidities of the negatively and positively charged lepton.In events with small values of φ * the pseudomass better approximates the dimuon mass.The pseudomass distributions for events with φ * < 0.05 are studied in intervals of φ and η of the µ ± candidate, with a further categorisation into candidates traversing the silicon strip or straw drift tube detectors downstream of the magnet.A maximum likelihood fit of the M ± distributions is performed for each of these detector regions.The signal shapes are described by the sum of a resonant Crystal-Ball [36] component and a nonresonant component represented by an exponential function.The means of the M ± Crystal-Ball functions are parameterised as M (1 ± A), where A and M are freely varying asymmetry and mass parameters, respectively.The resulting q/p corrections, which are given by A/p where p is the average muon momentum for a given interval in η and φ, are presented for the data and simulation for both polarity configurations in Fig. 2.
After the curvature corrections are applied to the data and simulation, the momenta of the simulated muons are smeared to match those in the data, as described below, according to where N (a, b) represents a random number sampled from a Gaussian distribution with mean a and width b.The σ MS and σ δ parameters correspond to the multiple scattering and curvature measurement contributions to the resolution, respectively.The smearing model includes six parameters in total.There are two momentum scale parameters α corresponding to the 2.2 < η < 4.4 region, which coincides with the selection of W boson candidates, and the η < 2.2 region.A single δ parameter, corresponding to a curvature bias, covers the region 2.2 < η < 4.4, while the value of δ is fixed to zero in the region η < 2.2.There are two σ δ parameters corresponding to the 2.2 < η < 4.4 and η < 2.2 regions, while a single σ MS parameter is found to adequately cover all η values.The empirical 1/cosh η dependence of the second term in Eq. 5 improves the modelling of the η dependence in the Z → µµ mass distribution.As a further correction to Eq. 5, the value of σ MS is increased by a factor of 1.5 in the region η > 3.3 since this improves the agreement between data and simulation in the η dependence of the quarkonia mass distributions.
The six smearing parameters are determined in a simultaneous fit of J/ψ → µµ, Υ (1S) → µµ and Z → µµ candidates in data and simulation.A total of 36 dimuon invariant mass distributions are used in the fit.First, there are three η regions covering η < 2.2, 2.2 < η < 3.3 and 3.3 < η < 4.4, which result in six categories that depend on the η regions of the two muons.The quarkonia mass distributions are only used in categories with both muons having η > 2.2.In the subset of the η categories with both muons in η > 2.2, the Z → µµ data are split into three intervals of the asymmetry between the

Parameter
Fit value momenta of the two muons, which provides a first order sensitivity to the δ parameters.
Finally, all categories are divided by magnet polarity.
As in previous studies of Z → µµ production with the LHCb experiment [37], the background under the Z → µµ peak is low enough to be neglected but the fit includes exponential functions for the background contributions under the quarkonia resonance peaks.The fractions and slopes of these exponential components vary freely in the fit.
The total χ 2 from the fit is 1862 for 2082 degrees of freedom.Table 1 shows the fit values of the six parameters in the smearing model.The δ values are close to zero as expected given the curvature corrections that have already been applied.Figure 3 shows the dimuon mass distributions for the Z, Υ (1S) and J/ψ samples after combining all categories with both muons in the 2.2 < η < 4.4 range.
The statistical uncertainties in the smearing parameters result in an uncertainty in m W of 3 MeV.The uncertainty in the world average of the Υ (1S) mass [7] leads to an uncertainty of 2 MeV.The uncertainty in the J/ψ mass is negligible compared to that in the Υ (1S) mass.The Υ (1S) and the Z masses have comparable relative uncertainties but the latter has a negligible effect given the limited size of the Z boson sample.The amount of material in the detector, which affects the modelling of energy losses, is varied by 10%, leading to an uncertainty of 3 MeV.A total uncertainty of 5 MeV is attributed to the shifts in m W corresponding to: an alternative form of the 1/cosh η factor in Eq. 5; and a variation in the η region over which σ MS is scaled.An uncertainty of 2 MeV is attributed to the modelling of the radiative tails in the Υ (1S) and J/ψ simulation, using the methods described in Ref. [38].The total uncertainty attributed to the modelling of the momentum scale and resolution is 7 MeV.

Efficiency corrections
Corrections to the simulation are required for the muon trigger, identification, tracking and isolation efficiencies.The efficiencies are measured using a combination of Z → µµ and, in the case of the trigger efficiency, Υ (1S) → µµ samples.Positively and negatively charged muons are analysed separately but the results are combined since any charge asymmetries are verified to have a negligible effect.
The trigger efficiency, which accounts for the hardware and software stages, is measured using a combination of Z → µµ and Υ (1S) → µµ events in which one so-called tag muon is required to match a positive decision in the hardware trigger and the first stage of the software-level trigger such that the other muon can be regarded as an unbiased probe of the trigger efficiency.Events are categorised as either matched or unmatched depending on whether the probe muon is matched to a positive trigger decision in the event data record.The Z → µµ sample is verified to be sufficiently pure that the efficiencies can be measured by simply counting the matched and unmatched events with invariant masses within ±15 GeV of the known Z boson mass [7].The efficiencies are determined in four uniform φ intervals and eight uniform η intervals in the range 2.2 < η < 4.4.There are two additional η intervals in the region η < 2.2 and one in the region η > 4.4.The Υ (1S) → µµ sample requires background subtraction by fitting the dimuon invariant mass distribution with a parametric model of the signal and background components.
Three p T intervals, in the range 7.0 < p T < 12.5 GeV, are used for the probe muons from Υ (1S) → µµ decays while for the Z → µµ candidates an adaptive algorithm is employed to determine the p T intervals.The ratios of the trigger efficiencies in data relative to those in the simulation are shown as a function of 1/p T of the muon in Fig. 4 for each of the intervals in η and φ.These are overlaid with a linear function of p T , from which correction weights for the simulated events are evaluated.The weights for the W boson model only rely on these functions but the weights for the Z boson model also require a parameterisation of the absolute efficiency in the simulation such that the efficiency can be correctly modelled for Z boson candidates with one or two muons matched to a trigger decision.The absolute efficiency is described by an error function that captures the p T threshold (roughly 6 GeV) of the hardware trigger [28].The muon identification efficiency is treated in a similar manner to the trigger efficiencies, using Z → µµ events.The resulting event weights, which are applied to the simulated events, are within a few per cent from unity.The tracking efficiency is determined as in previous measurements of W and Z boson production at LHCb [37] using Z → µµ candidates where the probe muons are reconstructed by combining hits from the muon stations and the large-area silicon-strip detector located upstream of the magnet [39].As neither of these detectors are used in the primary track reconstruction algorithms, the probes can be used to measure the tracking efficiency.Correction factors are evaluated using a similar approach to those of the muon identification efficiency, except that the corrections are assumed to be independent of p T .
The statistical uncertainties in the muon trigger, tracking and identification efficiency corrections are evaluated by rerunning the relevant steps of the analysis, up to and including the m W determination, with random fluctuations in the underlying efficiency values.The RMS of the resulting variations in the m W value is regarded as an uncertainty.A systematic uncertainty is attributed to the dependence of the results on the scheme for η and φ intervals.This includes restricting to a single interval in φ, reducing the number of η intervals (within 2.2 < η < 4.4) by a factor of two, varying the number of p T intervals between two and ten, and using the simulation rather than the data to control the adaptive algorithm.Further systematic uncertainties are attributed to variations in the isolation and p T requirements on the tag muons, the mass windows used to determine  the Z → µµ signal yields, and the functions of the p T dependence of the efficiency ratios.As the probe muons for the tracking efficiency are reconstructed using minimal tracking information, they have a significantly lower momentum resolution and so a dedicated momentum smearing is applied by default.A variation in the size of this smearing is included in the systematic uncertainty evaluation.The total uncertainty associated to the muon trigger, identification and tracking uncertainties is 6 MeV.
The efficiency of the isolation requirement is measured with Z → µµ events.The isolation variable receives contributions from pile-up, the underlying event and the recoil component of the hard process.The recoil projection for each muon is defined by where p µ T and p V T are the two-dimensional momentum vectors of the muon and the parent vector boson in the transverse plane.Figure 5 (left) shows the isolation requirement efficiency as a function of u in the Z → µµ data and simulation.The efficiency is around 80% at positive values of u where the underlying event contribution dominates.At negative values of u, corresponding to large recoil, the efficiency drops to around 70%.The full reconstruction of the Z boson in Z → µµ events allows the determination of corrections as a function of u defined at reconstruction level and consistently applied to W and Z boson events as a function of u defined at generator level, with this approach validated using the Z boson sample.A map of relative efficiencies between data and simulation is determined in intervals of u and η and is used to evaluate weights for the simulated events.Figure 5 (right) shows that the dependence of the isolation efficiency on 1/p T in data is accurately described by the simulation after the corrections.The statistical uncertainties in the isolation efficiency corrections are treated as a source of systematic uncertainty.This is combined in quadrature with a systematic uncertainty that accounts for variations in the u and η intervals and in a smoothing procedure applied to enhance the effective statistical precision of the correction map.The total uncertainty attributed to the isolation efficiency modelling is 4 MeV.
The modelling of the impact parameter and track fit χ 2 variables in simulation is improved in two stages, both of which make use of Z → µµ events.Initially, the values of the variables in the simulation are smeared and shifted to match the data.Subsequently, weights are applied to the simulated events to correct for small residual differences in the efficiencies of the selection requirements between data and simulation.In order for the impact parameter modelling to be reliably transported between Z and W boson events, the PVs are refitted with all signal muons removed.The three-dimensional impact parameter is then decomposed into its individual components and these are smeared according to a normal distribution in six intervals in η and seven intervals in φ.A similar procedure is used to improve the modelling of the track fit χ 2 distribution.These corrections are followed by smaller corrections applied to account for the efficiencies of the impact parameter and track χ 2 requirements.The efficiency weights are also determined with Z → µµ events and are typically within a few per mille from unity.Neither the statistical uncertainties nor reasonable variations in the η and φ interval schemes are found to have significant impact on the m W value.Therefore, no systematic uncertainty is considered.

QCD background model
A small background from in-flight decays of pions and kaons into muons is present in the sample of W → µν boson candidates.This background cannot be modelled with high enough accuracy using full detector simulation.It is therefore modelled using a sample of high-p T tracks, selected by dedicated triggers without muon identification requirements.The W boson selection requirements are applied to this sample but with the muon identification requirement inverted.The resulting sample is verified in simulation to be a pure sample of charged hadrons, composed of roughly 60% pions, 30% kaons and 10% protons, produced directly at the pp interaction vertex.In particular, the impact parameter requirements suppress the heavy flavour hadron content to a negligible level.
The probability of an unstable hadron of mass m, lifetime τ , and momentum p to decay within a detector of length d is Similar kinematic distributions are predicted for pions, kaons and protons in the simulation.Therefore, the in-flight decay background can be modelled by the data with weights of 1/p.The majority of the in-flight decays occur outside the magnetic field region and therefore have minimal influence on the measured momentum.The absolute normalisation is not needed because this background component is allowed to vary freely in the m W fit.
The weighted p T spectra for both charges are shown in Fig. 6 and are overlaid with best-fit functions of the form [40], where a and n are empirical parameters that are determined in the fits.In addition to giving a good fit to the data, this functional form is verified to describe the pion and kaon spectra in simulation.The fit functions are sampled to generate background candidates  for inclusion in the m W fit.A charge asymmetry of O(10%), favouring positively charged hadrons due to the pp initial state, is observed and is included in the sampling.
The uncertainty in the hadronic background model is dominated by systematic sources.Three different systematic uncertainties, which are combined in quadrature, are assessed.The first accounts for assuming that the data sample can be treated as containing a single hadron species.The second accounts for a small bias from the inverted muon identification requirements.The third accounts for a small dependence on the range of p T values used in the fits to Eq. 8.The combined systematic uncertainty is 2 MeV.

Modelling W and Z boson production
The emulation of different m W hypotheses is achieved by assigning event weights based on a relativistic Breit-Wigner function with a mass-dependent width.Further weight-based corrections are also applied to simulated W and Z boson events to improve upon the limited formal accuracy of Pythia.A weighting to model QCD effects is applied in the basis of Eq. 2. A further weighting to model QED effects is applied as a function of the logarithm of the relative energy difference between the dimuon system before and after QED final-state radiation.
Higher order electroweak corrections are not included in the model.Instead, an uncertainty of 5 MeV is attributed to these missing corrections using samples of events generated using POWHEGBoxV2 [41][42][43] with and without electroweak corrections.These events are interfaced with Pythia and the uncertainty is evaluated using the data challenge methods described in Sect.8.1.

Candidate QCD programs
Five software programs, or combinations of these, are evaluated as potential candidates for the weighting of the simulated W and Z boson events.
1. Pythia: Events are generated using Pythia version 8.235 [17] with several values of the intrinsic transverse momentum (k intr T ) of the initial state partons and α s , closely following the work of Ref. [19].The NNPDF23 lo as 0130 qed [44] PDFs are used in the event generation.The events are weighted using the methods described in Ref. [45] to the NNPDF31 lo as 0118 [46] and CT09MCS [47] PDFs.
2. POWHEGPythia: Events are generated using POWHEGBoxV2 [48] with the NNPDF31 nlo as 0118 PDFs and are subsequently showered with Pythia version 8.244 [17].The event generation with POWHEGBoxV2 is repeated with different values of α s .The default Monash [18] tune of Pythia is used but event samples are generated with different values of k intr T , and with the same value of α s as used in POWHEGBoxV2.This results in a grid of predictions with different α s and k intr T values.

POWHEGHerwig:
Events are generated equivalently to those from POWHEGPythia but substituting Pythia with Herwig [49] for the partonshower stage.
4. Herwig: These events are also equivalent to those of POWHEGPythia except that the hard process and the parton shower are both fully implemented in Herwig [49].

DYTurbo:
The cross-sections and angular coefficients are computed at O(α 2 s ) accuracy using DYTurbo [50] with the NNPDF31 nnlo as 0118 PDFs [46].Predictions for the unpolarised cross-section include resummation to next-to-next-to-leading logarithms and are produced with several values of the g parameter that controls nonperturbative effects.
Histograms of the unpolarised cross-section in Eq. 2 and the angular coefficients are produced for all combinations of programs and tuning parameters.These histograms, which are used to determine event weights, have intervals in the transverse momentum, rapidity and mass of the vector boson.

QCD weighting and transverse momentum model
The simulated samples described in Sect. 2 can be weighted in the full five-dimensional phase space of vector boson decays, according to Eq. 2, to provide predictions based on different models of QCD.For the unpolarised cross-section such weights are found by interpolating between the generated histograms described above.
A detailed measurement of the angular coefficients in pp → Z → µµ at √ s = 8 TeV was reported by ATLAS [51].Predictions based on parton showers are generally found to be unreliable in predicting the angular coefficients.However, the ATLAS data are reasonably well described by O(α 2 s ) predictions from DYNNLO [52], on which DYTurbo is based.An exception is the difference between A 0 and A 2 , for which O(α 2 s ) is effectively only leading order, but the present measurement has a negligible sensitivity to this particular detail.Hereafter DYTurbo is used in the modelling of the angular coefficients.
Since the prediction of each angular coefficient relies on separate numerator and denominator calculations, there are four independent renormalisation and factorisation scales that are varied to assess the uncertainty associated with missing higher orders in α s .In Ref. [53] it is argued that fully correlating the scale variations between the numerator and denominator, which leads to a large degree of cancellation, may result in inadequate uncertainty coverage.The present analysis therefore follows the recommendation of Ref. [53], which is to vary the four scales independently by factors of 1  2 and 2 with the constraint that all ratios that could be constructed from the four scales are between 1  2 and 2. This results in an envelope of 31 values of m W that sets the associated uncertainty.
Figure 7 compares the p Z T distribution in the data with the Pythia simulation weighted to the different unpolarised cross-section predictions before and after tuning them to the data.Table 2 lists the χ 2 and the preferred parameter values for the fits with each model.DYTurbo gives a reasonable prediction but overestimates the number of events with large p Z T even with tuning of the g parameter.A reasonable initial description is provided by Pythia, which is to be expected since it has already been tuned to p Z T data.The POWHEGPythia, POWHEGHerwig and Herwig predictions poorly describe the shape of the p Z T distribution with their default values of α s = 0.118.Their descriptions of the p Z T distribution are greatly improved when their α s and k intr T parameters are tuned.Of these programs, POWHEGPythia gives the most reliable description of the data, with a preferred α s value of around 0.125.Large values of α s are also favoured by other models and in other studies of the p Z T distribution [54]. 5Therefore, POWHEGPythia with freely varying α s and k intr T values is selected for the default fitting model.The systematic uncertainty in the description of the p Z T and p W T shapes is evaluated with alternative predictions from: Pythia with the CT09MCS and leading-order NNPDF31 PDFs; Herwig and POWHEGHerwig with the next-to-leading-order NNPDF31 PDFs.The envelope of shifts in m W obtained from using these alternative descriptions is found to be 11 MeV, providing the dominant contribution to the systematic uncertainty associated with the modelling of the vector boson transverse momentum.
All of the event generator predictions can be weighted at leading-order to emulate event generation based on different PDFs.As discussed in Ref. [45] this weighting is not completely valid for events generated at next-to-leading order.Since POWHEGPythia allows in situ computation of next-to-leading order PDF weights, it is possible to directly estimate the inaccuracy of the leading-order approximation.With five NNPDF31 nlo as 0118 replicas it is verified that the differences between the leading-order and next-to-leadingorder weighting approaches are smaller than 1 MeV in the m W fit.
Since it is computationally expensive to determine fully the PDF uncertainty in the DYTurbo angular coefficients, the variations in the next-to-leading order PDFs used in the POWHEGPythia model of the unpolarised cross-section are coherently propagated to the angular coefficients.The DYTurbo angular coefficients are shifted by the differences in values predicted by POWHEGPythia in the default (NNPDF31 nlo as 0118) PDF compared to the target PDF in the uncertainty assessment.
Separate measurements of m W based on the NNPDF3.1 [46], CT18 [55] and MSHT20 [56] PDF sets, each with their own PDF uncertainty estimate, are reported.However, since these three sets are based on almost the same data, the central result of this analysis is a simple arithmetic average of the three results, under the assumption that the three PDF uncertainties are fully correlated.

Angular scale factors
The uncertainties in the angular coefficients from DYTurbo would lead to an uncertainty of O (30) MeV in m W , with the dominant contribution attributed to the A 3 coefficient.The   can be understood by inspection of Eq. 2: an increase in A 3 enhances the cross-section for events with large sin ϑ and cos ϕ.The contribution to the muon p T from the W boson mass scales with sin ϑ while the contribution from the transverse momentum of the W boson scales with ± cos ϕ for W ± boson production.By allowing a single A 3 scaling factor, which is shared between the W + and W − processes, to vary freely in the m W fit the angular coefficient uncertainty is reduced by roughly a factor of three, to 10 MeV.Effectively the resulting model only depends on DYTurbo for the kinematic dependence of A 3 , while all other coefficients are fully modelled by DYTurbo.

Parametric correction at high transverse momentum
While POWHEGPythia is shown in Sect.7 to describe the p Z T distribution in the region below 30 GeV, it systematically underestimates the cross-section at higher p Z T .This is expected due to the missing matrix elements for the production of a weak boson and more than one jet.Figure 8 compares the p Z T distribution in the data with the model prediction

Ratio
Figure 9: Logarithm of the relative energy loss of the dilepton system due to final-state radiation for (left) W boson events and (right) Z boson events.An energy loss of below 10 −6 is considered unresolvable and is accounted for in the underflow bin to the left of the dashed vertical line.In the lower panel the ratio with respect to Pythia is shown.having set α s and k intr T to be close to the final fit values.For p Z T ≥ 40 GeV the model starts to underestimate the cross-section, reaching the ten per cent level at p Z T ∼ 100 GeV.In the lower panel of Fig. 8 the data to prediction ratio is overlaid with a function of the form Since the universality of this correction between W and Z boson processes is not well controlled, an uncertainty of 100% of this correction is included as an additional systematic uncertainty associated with the vector boson p T model.This contributes an uncertainty in the m W value smaller than 1 MeV.

QED weighting
The effect of the QED final-state radiation is largely characterised by the energy difference between the final-state lepton system before and after radiation.The logarithm of this energy difference as described by Pythia, Herwig, and an alternative configuration of Pythia with the final-state radiation modelled by Photos [57] is shown in Fig. 9. Event-by-event weights are evaluated for the Herwig and Photos models relative to Pythia and applied to the simulated samples used in the analysis.The default model uses the arithmetic average of the Herwig, Photos and Pythia weights, where the Pythia weights are equal to unity.The systematic uncertainty, amounting to 7 MeV, is taken from the envelope of m W values corresponding to each of the three models taken individually.

W boson mass fit
The m W value is determined through a simultaneous fit of the q/p T distribution of W boson candidates and the φ * distribution of Z boson candidates.The φ * variable is preferred over p Z T because it is less susceptible to several of the modelling uncertainties, while still being sensitive to the parameters affecting the predicted p Z T and p W T distributions.The φ * distribution extends to φ * = 0.5 while the q/p T distribution includes two fit regions covering 28 < p T < 52 GeV.Projections of the q/p T distribution cover a wider interval with p T > 24 GeV that includes regions outside the fit.The model of both distributions is based on simulated event samples with event-by-event weights.The fit minimises the sum of two negative log-likelihood terms, associated with the q/p T and φ * distributions, that are computed using the Beeston-Barlow-Lite prescription [58], which accounts for the finite size of the simulated samples.The sum of these terms multiplied by a factor of two is denoted as the χ 2 .
The φ * distribution is modelled including background contributions from Z → τ τ and top quarks.The model of the q/p T distribution includes the dominant W → µν signal component and several background sources.The largest background, with a fraction of around 7 × 10 −2 , is attributed to Z → µµ, which is simulated with true invariant masses above 20 GeV.The W → τ ν and hadronic background components each contribute at the O(10 −2 ) level.A combination of rarer background sources including Z → τ τ , top quarks, vector boson pairs, and heavy flavour hadrons gives a total contribution below 10 −2 .
The fractions of the W + and W − signal components and the hadronic background are allowed to vary freely.The W → τ ν component is constrained, using the known τ → µν ν branching fraction [7], relative to the W → µν component.All other component fractions are fixed relative to the observed number of Z boson candidates in the φ * distribution using the fiducial cross-sections for the corresponding processes relative to that of Z boson production.The fiducial selections for W and Z boson processes are the same as used in measurements of the corresponding cross-sections [59].For all other processes the fiducial regions correspond to the requirement of a single muon in the region p T > 20 GeV and 2 < η < 4.5.The measured cross-section for Z → µµ is used [37].The cross-sections for the rare background processes are determined using Powheg with the next-to-leading-order NNPDF3.1 PDF sets.
The shapes of background components arising from the decay of electroweak bosons are determined from the same models used to describe the signal component.Systematic variations of the model used to describe W boson production therefore also simultaneously provide systematic variations associated with the shapes of background contributions from the decay of electroweak bosons.The uncertainty in m W from varying the predicted cross-sections for the rarer background within their uncertainties is negligible.
The default physics model is based on the simulated samples set out in Sect.2, fully weighted using a combination of POWHEGPythia and DYTurbo for the QCD description and a combination of Pythia, Photos and Herwig for the QED description.The fit is configured to determine the following parameters:

Data challenge tests
In Sect.7 it is concluded that POWHEGPythia describes the p Z T distribution, in the p Z T ≤ 30 GeV region, better than the other candidate models.It is important to demonstrate that the fit can reliably determine m W if W boson production is better described by one or more of the other models.Several pseudodata samples are prepared in which the underlying Pythia events, without detector simulation, are weighted to match the default DYTurbo and POWHEGPythia model but with the p V T distribution modified to match an alternative model.The m W fit is configured with a simplified model, without background components, using a statistically independent sample of the same Pythia events without detector simulation.Figure 10 shows the resulting q/p T and φ * distributions of these pseudodata samples.Variations of up to five per cent are seen in the shape of the φ * distribution.Within the fit regions in q/p T variations of several per cent can be seen, while the variations exceed O(10 −1 ) in the high-p T control region.However, the fit model is able to absorb these differences in the α s and k intr T nuisance parameters with variations in the preferred m W value of no more than 10 MeV.Table 3 lists the results of the fits to these pseudodata samples.The observed variation in m W is consistent with the uncertainty due to modelling the vector boson transverse momentum distribution in the fit to LHCb data, as discussed in Sect.7.

Fit results
The fit to the data, with the NNPDF31 nlo as 0118 PDF set, returns a total χ 2 of 105 for 102 degrees of freedom.Figure 11 compares the q/p T and φ * distributions from the data with the fit model overlaid.The model is in good agreement with the data within the fit ranges but it underestimates the high-p T control region of the q/p T distribution by up to ten per cent.This underestimation is within the band of modelling uncertainty, which is dominated by the high-p V T parametric correction in that region.The values of the eight  The DYTurbo pseudodata correspond to α s = 0.118 and g = 1 GeV 2 .The contributions to the total χ 2 from the q/p T and φ * distributions are denoted χ 2 W and χ 2 Z , respectively.The shift in the m W value with respect to the POWHEGPythia pseudodata is denoted δm W .The uncertainties quoted are statistical.parameters determined from the fit are listed in Table 4.The α s value for the W boson events is roughly 0.002 higher than for the Z boson events.If the fit is configured with a shared α s value for the W and Z boson events the value of m W changes by +39 MeV but the χ 2 is increased by more than 20 units, which strongly favours the configuration with independent α s values.Furthermore, similar variations between the α s values for W and Z boson events are found in the data challenge tests, as shown in Table 3.The A 3 scaling factor is statistically consistent with unity, which suggests that the O(α 2 s ) predictions from DYTurbo, with the central scale choices, are compatible with the data.
Figure 12 (left) shows the projection of the q/p T distribution in the Z boson sample, where the final state muon is only included if it satisfies the W boson selection requirements.The model is in good agreement with the data.Figure 12 (right) shows that the Z boson rapidity distribution is well described by the model.The uncertainties are evaluated according to the specific methods for the three groups.The NNPDF3.1 uncertainty is evaluated as the RMS of m W values according to 100 replicas, whereas the other two sets use fixed numbers of eigenvectors.The CT18 uncertainty is corrected from 90% confidence level (CL) to 68% CL to be consistent with all other uncertainties in this analysis.For each PDF set, the uncertainty from the replica variations is added in quadrature to the uncertainty from variations in the α s used in the PDF fits.Values of 0.116 < α s < 0.120 are considered, and the uncertainty in m W is taken as half of the absolute difference between the corresponding shifts in m W [60]. 6 The PDF uncertainty on the arithmetic average of the three results is taken as the arithmetic average of the three uncertainties, in accordance with the assumption that the uncertainties are fully correlated.Table 6 lists each contribution to the systematic uncertainty in the final result, after averaging results based on the three PDF sets.The systematic uncertainty is split into three orthogonal components that are combined in quadrature.The uncertainty due to the description of the parton distribution functions is 9 MeV.The remaining theory uncertainty in the modelling of W and Z boson production is 17 MeV, as described in Sect.7, with the largest contribution arising from variations of the transverse momentum model.The experimental uncertainty is 10 MeV, with the different contributions discussed in Sects.4, 5, and 6.
Independently of the systematic uncertainty evaluation, several cross-checks of the measurement are performed.
• Consistency of orthogonal subsets: The data and simulation are split into orthogonal subsets by magnet polarity, the product of the muon charge and polarity, and the φ and η of the muon in the W boson selection.These results are reported in Table 7. Considering the statistical uncertainties only, all differences are within, or just outside, two standard deviations, which was predefined as a criterion for this test.
• Fit range: The minimum and maximum p T of the fit range in the q/p T distribution • Fit model freedom: The choice of parameters that are determined in the fit is varied and the results are reported in Table 9.The default fit determines one α s parameter for the Z processes and a second that is shared between W + and W − processes.With three α s parameters there is only a small change in m W and the fit quality.The default fit determines a single floating k intr T parameter that is shared among all three processes.Neither the m W value nor the χ 2 are strongly affected by allowing two (with one shared between the W + and W − processes) or three k intr T parameters to vary freely.If the A 3 scaling factor is fixed to unity the value of m W shifts by 7 MeV and the χ 2 increases by a few units.In summary the m W fit seems to be rather insensitive to all of these variations, except that the data strongly prefer independent POWHEGPythia tunes for the W and Z boson production processes.
• Use of NNLO PDF sets: The PDF set used for the analysis is varied from NNPDF31 nlo as 0118 to NNPDF31 nnlo as 0118.The shift in m W is 1 MeV.
• Separate m W values for W + and W − bosons: an additional parameter is included in the fit, allowing for separate values of m W for W + and W − bosons.This mass difference is found to be consistent with zero within one standard deviation.
• W -like measurement of the Z boson mass: the same methods are applied to the Z boson sample alone, to perform a W -like measurement of the Z boson mass.The values measured with positive and negatively charged muons agree within one standard deviation and their average is consistent with the PDG average [7] within one standard deviation.[61], DELPHI [62], L3 [63], OPAL [64], CDF [10], D0 [11] and ATLAS [12] experiments.The current prediction of m W from the global electroweak fit is also included.

Summary and Conclusion
This paper reports the first measurement of m W with the LHCb experiment.A data sample of pp collisions at √ s = 13 TeV corresponding to an integrated luminosity of 1.7 fb −1 is analysed.The measurement is based on the shape of the p T distribution of muons from W boson decays.A simultaneous fit of the q/p T distribution of W boson decay candidates and of the φ * distribution of Z boson decay candidates is verified to reliably determine m W .This method has reduced sensitivity to the uncertainties in modelling the W boson transverse momentum distribution compared to previous determinations of m W at hadron colliders.The following results are obtained m W = 80362 ± 23 stat ± 10 exp ± 17 theory ± 9 PDF MeV, m W = 80350 ± 23 stat ± 10 exp ± 17 theory ± 12 PDF MeV, m W = 80351 ± 23 stat ± 10 exp ± 17 theory ± 7 PDF MeV, with the NNPDF3.1,CT18 and MSHT20 PDF sets, respectively.The first uncertainty is statistical, the second is due to experimental systematic uncertainties, and the third and fourth are due to uncertainties in the theoretical modelling and the description of the PDFs, respectively.Treating the three PDF sets equally results in the following arithmetic average m W = 80354 ± 23 stat ± 10 exp ± 17 theory ± 9 PDF MeV.
This result agrees with the current PDG average of direct measurements [7] and the indirect prediction from the global EW fit [6], and is compared to previous measurements in Fig. 13.This measurement also serves as a first proof-of-principle of a measurement of m W with the LHCb experiment.In Ref. [65] it was demonstrated that the PDF uncertainty in a measurement of m W by LHCb can be strongly reduced by using in situ constraints and by fitting the doubly differential distribution of p T and η, similar to the measurement by the CMS Collaboration [66], instead of the singly differential p T distribution.An approximately three times larger data sample is already available for analysis but particular attention should be paid to reducing the dominant source of systematic uncertainty, which is the modelling of W boson production.

Figure 1 :
Figure 1: Muon q/p T distribution in simulated W → µν events with variations in (left) m W and (right) the A 3 coefficient.The dashed vertical lines indicate the two p T ranges that are included in the m W fit.

Figure 2 :
Figure 2: Curvature corrections as a function of the detector region index (depends on η, φ and tracking detector, as described in the text) for (left) data and (right) simulation.The corrections are shown for both polarity configurations.The periodic pattern corresponds to a dependence on pseudorapidity that repeats in the intervals of the azimuthal angle.

Figure 3 :
Figure 3: Dimuon mass distributions for selected J/ψ, Υ (1S) and Z boson candidates.All categories with both muons in the 2.2 < η < 4.4 region are combined.The data are compared with the fit model.The red histogram delineates the model before the application of the smearing.

Figure 4 :
Figure4: Trigger efficiency ratios in data relative to simulation in intervals of η for a single interval of φ.The points are presented as a function of 1/p T and are overlaid with best-fit linear functions of p T .Points from each η interval are separated for readability by an offset of half the integer index i for that interval.The uncertainties on most of the points are too small for the error bars to be visible.

Figure 5 :
Figure 5: Isolation efficiencies as a function of the observables (left) u and (right) 1/p T of the muon.The grey histograms indicate the, arbitrarily normalised, shapes of each distribution in simulated W boson events.In the lower panels the ratios of the isolation efficiency with respect to the uncorrected simulation are shown.

Figure 6 :
Figure 6: Weighted p T spectra of the samples of (left) positively and (right) negatively charged hadron candidates.The fit results are overlaid.

Figure 7 :
Figure 7: Distributions of p Z T (left) before and (right) after the fit for the different candidate models of the unpolarised cross-sections.The fit only considers the region p Z T < 30 GeV, indicated by the dashed vertical line.In the lower panels the ratios with respect to the POWHEGPythia model are shown.

Figure 8 :
Figure 8: Distribution of p Z T compared to the POWHEGPythia model prior to the parametric correction, which is delineated by the red line in the lower panel.

1 .
the value of m W , 2. the fraction of W + signal, 3. the fraction of W − signal, 4. the fraction of QCD background, 5. the value of α s for the Z boson processes (α Z s ), 6. an independent α s value that is shared for the W + and W − signals (α W s ), 7. a shared k intr T value for all W and Z boson processes, 8. and an A 3 scale factor that is shared by the W + and W − signals.

Figure 10 :
Figure 10: Projections of the (left) q/p T and (right) φ * distributions for the challenge datasets.The four dashed vertical lines indicate the two fit regions in the q/p T distribution.

Figure 11 :
Figure 11: Distributions of (left) q/p T and (right) φ * compared to the model after the m W fit.

Figure 12 :
Figure12: Projections of the (left) q/p T and (right) rapidity distributions for the Z boson selection.A final state muon is only included in the q/p T distribution if it satisfies the W boson selection requirements.

Table 1 :
Parameters in the momentum smearing model where the uncertainties quoted are statistical.

Table 2 :
Results of fits of different models to the p Z T distribution.The uncertainties quoted are statistical, and the χ 2 comparison of the different models to the data is evaluated considering only statistical uncertainties.The right-hand column lists the fit values of the k intr T parameter or, for DYTurbo, the analogous g parameter.The fit with DYTurbo has one more degree of freedom than the fits with the other models since only one tuning parameter (g) is used for DYTurbo.

Table 4 :
Values of the parameters determined in the m W fit with the NNPDF31 nlo as 0118 PDF set.The uncertainties quoted are statistical.

Table 5 :
Uncertainties for the NNPDF3.1,CT18 and MSHT20 sets.The contributions from the PDF uncertainty with fixed α s and from the α s variation are quoted separately as is their sum in quadrature, which defines the total uncertainty for each PDF set.

Table 5
lists the PDF uncertainties evaluated for fits based on the NNPDF3.1,CT18 and MSHT20 PDF sets.The m W values agree within an envelope of 12 MeV, which supports the choice to report an arithmetic average of the three.With respect to the central value obtained with the NNPDF31 PDFs, the m W values obtained with the CT18 and MSHT20 PDF sets differ by −12 MeV and −11 MeV, respectively.

Table 6 :
Contributions to the systematic uncertainty in m W . Negligible contributions below 1 MeV are not listed.theirdefaultvalues of 28 GeV and 52 GeV, respectively.The results are reported in Table8.Considering the variations in the statistical uncertainties in m W this test shows that the fit results are stable with respect to variations in the fit range.

Table 7 :
Fit results where the data and simulation samples are split into two orthogonal subsets.For a given split, the first row is defined as the reference with respect to which the difference in m W , denoted by δm W , is defined.The uncertainties quoted on δm W are statistical.

Table 8 :
Fit results with variations in the fit range around the default p min T = 28 GeV and p min T = 52 GeV.The second column lists the χ 2 values, the third column lists the shifts in m W with respect to the default fit and the third column lists the statistical uncertainties in m W .

Table 9 :
Fit results with variations in which physics parameters are varying freely.
Universidade Federal do Triângulo Mineiro (UFTM), Uberaba-MG, Brazil b Hangzhou Institute for Advanced Study, UCAS, Hangzhou, China c Università di Bari, Bari, Italy d Università di Bologna, Bologna, Italy e Università di Cagliari, Cagliari, Italy f Università di Ferrara, Ferrara, Italy g Università di Firenze, Firenze, Italy h Università di Genova, Genova, Italy i Università degli Studi di Milano, Milano, Italy j Università di Milano Bicocca, Milano, Italy k Università di Modena e Reggio Emilia, Modena, Italy l Università di Padova, Padova, Italy m Scuola Normale Superiore, Pisa, Italy n Università di Pisa, Pisa, Italy o Università della Basilicata, Potenza, Italy p Università di Roma Tor Vergata, Roma, Italy q Università di Siena, Siena, Italy r Università di Urbino, Urbino, Italy s MSU -Iligan Institute of Technology (MSU-IIT), Iligan, Philippines t AGH -University of Science and Technology, Faculty of Computer Science, Electronics and Telecommunications, Kraków, Poland u P.N.Lebedev Physical Institute, Russian Academy of Science (LPI RAS), Moscow, Russia v Novosibirsk State University, Novosibirsk, Russia w Department of Physics and Astronomy, Uppsala University, Uppsala, Sweden a x Hanoi University of Science, Hanoi, Vietnam