Luminosity Spectrum Reconstruction at Linear Colliders

A good knowledge of the luminosity spectrum is mandatory for many measurements at future e+e- colliders. As the beam-parameters determining the luminosity spectrum cannot be measured precisely, the luminosity spectrum has to be measured through a gauge process with the detector. The measured distributions, used to reconstruct the spectrum, depend on Initial State Radiation, cross-section, and Final State Radiation. To extract the basic luminosity spectrum, a parametric model of the luminosity spectrum is created, in this case the spectrum at the 3 TeV Compact Linear Collider (CLIC). The model is used within a reweighting technique to extract the luminosity spectrum from measured Bhabha event observables, taking all relevant effects into account. The centre-of-mass energy spectrum is reconstructed within 5% over the full validity range of the model. The reconstructed spectrum does not result in a significant bias or systematic uncertainty in the exemplary physics benchmark process of smuon pair production.


Introduction
Small, nanometre-sized beams are necessary to reach the required luminosity at future linear colliders. Together with the high energy, the small beams cause large electromagnetic fields during the bunch crossing. These intense fields at the interaction point squeeze the beams. This so-called pinch effect increases the instantaneous luminosity. However, the deflection of the particles also leads to the emission of Beamstrahlung photons -which reduce the nominal energy of colliding particles -and produces collisions below the nominal centre-of-mass energy [1,2,3,4]. The resulting a e-mail: stephane.poss@cern.ch b e-mail: andre.sailer@cern.ch spectrum of centre-of-mass energies is traditionally called the luminosity spectrum 1 [4,5,6,7].
The knowledge of the shape of this luminosity spectrum is mandatory for the precision measurements in which a cross-section has to be known. While the cross-section depends on the centre-of-mass energy, the observables measured in the lab frame also depend on the difference in energy of the colliding electrons 2 , which determines the Lorentz boost of the system.
Unlike the electron structure functions (i.e, Initial State Radiation (ISR)) -which can be calculated precisely -the beam-beam forces, and therefore the Beamstrahlung, highly depend on the geometry of the colliding bunches. The actual beam-beam interaction taking place at the interaction point cannot be precisely simulated, because the geometry of the bunches cannot be measured. Therefore, the luminosity spectrum at the interaction point has to be measured using a physics channel with well known properties, e.g., Bhabha scattering.
The observables measured in the events are affected by detector resolutions. The distributions used for the reconstruction of the luminosity spectrum are also dependent on the cross-section of the process, and Initial and Final State Radiation (FSR). All effects have to be taken into account for the reconstruction of the luminosity spectrum.
It was pointed out by Frary and Miller [5] that a precise reconstruction of the peak of the luminosity spectrum, necessary for a top-quark threshold scan, can only be achieved with a measurement of the angles of the outgoing electrons from Bhabha scattering. The angles of the two particles are the most precisely measurable observable [5]. The angles of the outgoing electrons -or rather the acollinearity between the two particles -are sufficient to extract a relative centre-of-mass energy, which gives access to the luminosity spectrum.
Toomi et al. [6] showed that the reconstruction of a parameterised luminosity spectrum is possible using a template fit. Their parameterisation only used three parameters to describe the effective centre-of-mass energy spectrum. However, as the boost of the initial system and correlation between the energies of the two particles cannot be neglected [8], a description of the energies of the pairs of colliding particles is necessary. Correlations exist between the two particle energies because the probability to emit beamstrahlung depends on the distance travelled in the field of the opposite bunch, and the field strength depends on the position inside the bunch. As two particles can only collide, when they are in similar position in their respective bunches the energy between two particles is correlated.
The relative centre-of-mass energy, that is reconstructed from the acollinearity, is equal to unity for back-to-back particles, and always smaller than unity for larger acollinearity, regardless whether one of the particles has a higher or lower energy than nominal. Therefore, Shibata et al. [7] proposed to calculate the distribution of the four-vectors of the Bhabha electrons and extract the luminosity spectrum with the iterative Expectation-Maximisation algorithm. They have, however, considered neither detector resolutions, nor Initial and Final State Radiation. For a full description of the outgoing Bhabha electrons, the luminosity spectrum would have to be weighted with the Bhabha cross-section and convoluted with the detector resolutions, which would require a huge computational effort, when using their method.
A reconstruction of the energy of the particle pairs was done for the 500 GeV ILC [9]. The acollinearity and the energies of the electrons measured in the calorimeter were used in a reweighting fit to reconstruct the luminosity spectrum. The parameterisation -necessary for the reweighting fit -accounted for the correlation between the two beams and the beam-energy spread.
This paper follows the approach of the 500 GeV ILC study [9], extends it, and applies it to the luminosity spectrum of the 3 TeV CLIC [10], which is the most challenging luminosity spectrum. This paper is structured as follows: in Section 2 the basic and cross-section scaled luminosity spectrum are defined. The Bhabha scattering and observables used for the reconstruction are also introduced. In Section 3 the model of the luminosity spectrum, required to perform a reweighting fit, is constructed. The reweighting technique is explained in Section 4, and in Section 5 it is applied to first validate the model against the luminosity spectrum at the 3 TeV CLIC; then all the relevant effects leading to the measured observables are included, and the luminosity spectrum is reconstructed from these distributions. In Section 6 the impact of the reconstructed luminosity spectrum on the measurement of the masses of supersymmetric particles in a CLIC benchmark process is estimated. The paper closes with a summary, conclusions, and outlook in Section 7.
2 Luminosity Spectrum, Bhabha Scattering, and the Measurement The nominal centre-of-mass energy √ s nom of a collider with two beams with the nominal beam energy E Beam is given by √ s nom = 2E Beam . If the two interacting particles carry only a fraction of the nominal beam energy x 1,2 = E 1,2 /E Beam , the effective centre-of-mass energy is The basic luminosity spectrum L x describes either the distribution of the fraction of centre-of-mass energies x = √ s / √ s nom . or the distribution of the fraction of energies of colliding particles L x 1 , x 2 prior to hard collisions and prior to Initial State Radiation. The two functions are connected via the integral along the lines of constant centre-ofmass energies, given by Equation (1). Therefore, The luminosity spectrum affects all centre-of-mass energy dependent observables. For example, the luminosity spectrum has to be used to predict the inclusive (i.e., observed) cross-section σ Machine Eff at the machine. The principle is the same as for the parton density functions at hadron machines, except that the luminosity spectrum depends on the machine and not only on the colliding particles. To calculate the effective cross-section the differential cross-section is weighted with the luminosity spectrum, either with the one-dimensional luminosity spectrum or for the two-dimensional luminosity spectrum Bhabha scattering is the process of choice for luminosity measurements. It can be calculated with high precision and has a large cross-section. To first order, the differential Bhabha cross-section is [11] dσ Bhabha dθ = 2πα 2 s sin θ sin 4 (θ /2) , where α is the fine-structure constant and θ the polar scattering angle. Because cross-sections σ x √ s depend on the centreof-mass energy, any process used to reconstruct the basic luminosity spectrum will inherently contain a scaled luminosity spectrum This means that it is not enough to reconstruct the observed centre-of-mass energy spectrum, the luminosity spectrum has to be extracted from the observed spectrum.
The observed centre-of-mass energy is further affected by Initial State Radiation. It is impossible to distinguish between energy loss from Initial State Radiation and Beamstrahlung on an event-by-event basis. Initial State Radiation and Beamstrahlung have to be disentangled statistically.
Finally, the scattered particles are recorded in the detector, where their properties are reconstructed within the limits of the resolution of the respective sub-detectors.

The Basic Luminosity Spectrum
The luminosity spectrum distribution can be seen as a convolution of the beam-energy spread, which is inherent to the accelerator, and the Beamstrahlung due to the Beam-Beam effects. Figure 1 shows the beam-energy spread of the 3 TeV CLIC machine. It is obtained from a simulation of the main linear accelerator and the beam delivery system [12]. The energy of a particle depends on its longitudinal position in the bunch (Figure 1a). Due to intra-bunch wakefields, particles in the front of the bunch gain more energy from the RF cavities than particles in the back of the bunch [13]. This leads to the two distinct peaks near the minimal and the maximal value of the beam-energy spread ( Figure 1b). The energy spread is not following a Gaussian distribution. The dependence of the particle energy on the longitudinal position also leads to larger correlations between the particle energies. The precise shape of the beam-energy spread depends on the RF-phase and the bunch length [13]. To avoid a loss in the luminosity, these parameters are not allowed to vary freely and have to be precisely controlled [10]. Therefore, it can be assumed that a limited knowledge of the shape of beam-energy spread is available.
The distribution of particles is used as the input to the beam-beam simulation. The simulation of the beam-beam effects is done with GUINEAPIG [4]. During the bunch crossing the intense electromagnetic fields -due to the opposing bunches -deflect the beam particles and cause Beamstrahlung. Figure 2 shows the full range and the region around the maximal energy of the two-dimensional luminosity spectrum. The square region in the distribution of the two energies is due to the beam-energy spread (see Figure 1b). Events with x 1 < 0.995 or x 2 < 0.995 were significantly affected by the Beamstrahlung. Figure 3a shows the basic luminosity spectrum with respect to the effective centre-of-mass energy √ s . The spectrum possesses a peak around the nominal centre-of-mass energy and a long tail down to less than 5% of the nominal centre-of-mass energy. Figure 3b shows the peak of the luminosity spectrum as it is produced by GUINEAPIG. Because the beam-energy spread is not Normally distributed, the centre-of-mass energy peak is not Gaussian either. Figure 3b also shows a spectrum obtained by randomly pairing the energies of two particles, i.e., removing the correlation between the energies of the two beams. There is a clear difference between the two cases. If the correlation between the particle energies is not taken into account, the luminosity spectrum cannot be reconstructed properly.
To describe the beam-energy spread -and anchor the luminosity spectrum -the absolute energy of the beam has to be known. The average beam energy can be measured on a level of 0.04% [10] with a dipole and beam position monitor in the beam delivery system of the accelerator. If the distribution itself can be measured as well is still under study.

Cross-Section-Scaled Luminosity Spectrum
The observed events are distributed according to the scaled luminosity spectrum (Equation (6)), thus the events obtained from GUINEAPIG have to be either weighted with very large weights, or sampled with an accept-reject method [14] to obtain events with constant weights. Because large eventweights are undesirable, the accept-reject method is chosen.
The 3 TeV CLIC luminosity spectrum extends over more than three orders of magnitude of the Bhabha cross-section, meaning the accept-reject method is very inefficient. If a very large number of events for the basic luminosity spectrum were available, the scaled luminosity spectrum could be directly sampled from them. To avoid storing the large number of basic events the accept-reject method is directly added in GUINEAPIG.
The differential cross-section of the Bhabha scattering has to be known for the accept-reject method. Instead of using Equation (5) to calculate the Bhabha cross-section, it is estimated with BHWIDE [15]. BHWIDE includes higherorder effects and Initial State Radiation. Only events with the electron and positron polar angle inside the tracking acceptance (7 • < θ < 173 • ) are accepted. The cross-section is estimated at precise centre-of-mass energies from 10 GeV to  3000 GeV without any luminosity spectrum. Figure 4 shows the cross-section as given by BHWIDE. Figure 4 shows the basic luminosity spectrum obtained with GUINEAPIG, the bin-wise multiplication of the luminosity spectrum with the cross-section, and the scaled luminosity spectrum from GUINEAPIG with the cross-section used in the accept-reject method. The last two curves are nearly identical showing that the modified GUINEAPIG produces a properly scaled luminosity spectrum with equally weighted events.

Generation of Bhabha Events
The different luminosity spectra are used in the Bhabha generator to create the events which are observed in the detector. The Bhabha events are generated with BHWIDE, where the energies of the initial electron and positron can be defined on an event-by-event basis, as implemented by Rimbault et al. [16]. The polar angle θ of the final state electrons must be 7 • < θ < 173 • to ensure they will be observable in the tracker. BHWIDE produces also Initial and Final State Radiation photons and accounts for their effects during the Bhabha scattering. Bhabha cross-section GP Basic cross-section × GP GP Scaled Fig. 4 Bhabha cross-section from BHWIDE for electrons with a polar angle 7 • < θ < 173 • , and the basic luminosity spectrum from GUINEAPIG (GP Basic), the luminosity spectrum scaled by bin-wise multiplication with a normalised Bhabha cross-section (GP × cross-section), and the luminosity spectrum scaled with an accept-reject method in GUINEAPIG (GP Scaled). The lines of 'GP × cross-section' and 'GP Scaled' are overlapping.
The cross-section -including the luminosity spectrum -for events with a centre-of-mass energy above 1.5 TeV is 11 pb, which results in more than 1 million events for an integrated luminosity of 100 fb −1 . The expectation for the 3 TeV CLIC is an integrated luminosity of 500 fb −1 per year [17].

Observables and Detector Resolutions
Three observables, which can be extracted from the final state electrons, are used for the reconstruction of the luminosity spectrum: the relative centre-of-mass energy calculated from the polar angles of the outgoing electron and positron √ s acol / √ s nom , the energy of the electron E 1 , and the energy of the positron E 2 . The relative centre-of-mass energy reconstructed from the acollinearity of the final state electrons is [8,9] √ s acol √ s nom = sin(θ 1 ) + sin(θ 2 ) + sin(θ 1 + θ 2 ) sin(θ 1 ) + sin(θ 2 ) − sin(θ 1 + θ 2 ) , where θ 1 is the polar angle of the electron and θ 2 that of the positron 3 . √ s acol is equal to the effective centre-of-mass energy √ s , if only one of the particles radiated photons -Beamstrahlung or Initial State Radiation. If both the electron and the positron radiated photons, the reconstructed centreof-mass energy √ s acol will be larger than √ s . The GEANT4 simulation of tens of millions of electrons is too time-consuming. To include the detector effects, resolution functions of the energy and angles have been obtained from fully simulated and reconstructed Bhabha events using the CLIC ILD CDR detector model [18]. The dominating beam-induced background, the γγ → hadron events [19], was accounted for.
The rate of electrons produced in Bhabha scattering falls with an increasing polar angle θ (cf. Equation (5)) and the events will be predominantly at small polar angles. Because the magnetic field is nearly collinear to those tracks, their curvature does not allow for an accurate measurement of the momentum. Therefore, the energy is reconstructed using only the calorimeter information. The tracking information is used to measure the angles. The energy resolution is shown in Figure 5a. It is modelled in the analysis with obtained from the reconstruction with overlaid CLIC 3 TeV γγ → hadron background. The angular resolution needed for the computation of the relative centre-of-mass energy depends on the energy, shown in Figure 5b. The angular resolution is better than 20 µrad for particle energies above 200 GeV. The distributions of the particle energies and the relative centre-of-mass energy are shown in Figure 6. Figures 6a and  6b show the distributions before and Figure 6c and 6d after the application of the resolution effects via four-vector smearing. The relative centre-of-mass energy is hardly affected by the resolution, due to the high angular resolution of the tracking detectors. The energy of the particles is much more affected by the detector resolution.

Modelling the Luminosity Spectrum
For the reconstruction of the basic luminosity spectrum with the reweighting fit, a model or parameterisation of the luminosity spectrum is needed. If the beam energy were not affected by beam-energy spread or Beamstrahlung, it could simply be described by a Dirac delta-distribution with the random variates for this function x E Beam = 1. The nominal energy is modified by the contributions from beamenergy spread ∆ x Spread and Beamstrahlung ∆ x Strahl which change the beam energy away from its initial value. The functions describing the two contributions will therefore describe the difference to the nominal value, and the final random variate will be  Thus, the functions for the beam-energy spread and Beamstrahlung should describe the energy change of particles due to the respective effect.
The luminosity spectrum model will be built by describing the two particle energies. Using x 1 and x 2 as defined in Section 2, the simplest description of a two-dimensional However, a purely factorising Ansatz is insufficient to describe the correlation between the particle energies. Therefore, the two-dimensional energy distribution is divided into four regions (as shown in Figure 7): one region where neither particle radiated Beamstrahlung (called the 'Peak'); two regions where one or the other particle radiated Beamstrahlung (called the 'Arms'); and one region where both particles radiated Beamstrahlung (called the 'Body'). This separation is only determined by whether a particle produced Beamstrahlung or not and ignores the beam-energy spread for the moment. The result of this division is a piecewise function For each region, the resulting particle energies are described by a product of the functions for the two particles Region (x 2 ), and the individual functions are constructed from convolutions of the beam-energy spread and Beamstrahlung functions, depending on the region.

Parameterisation of the Beam-Energy Spread
The function f Peak to describe the behaviour of the beamenergy spread ( Figure 8a) has to rise very steeply at the two extremities. A hyperbolic cosine, parabola, or higher order polynomials -with a reasonable number of parameters -do not describe this energy distribution well. A beta- is used to describe the beam-energy spread. As the beta-distribution is limited between 0 and 1, a variable transform is used to describe the beam-energy spread between the two endpoints x min and x max near the maximal values at the beginning and end of the distribution. Here, the variable x is the relative difference between a particle's energy E Particle and the nominal beam energy E Beam , which corresponds to ∆ x Spread from Equation (10). To also describe the particles with energies below the x min and above x max , the beta-distribution is convoluted with a Gaussian distribution g x with a mean µ = 0 and a width σ . The Beam-Energy Spread (BES) function is where ⊗ is the convolution operator Due to Fubini's theorem the convolution of two probability density functions always results in another probability density function [20].
The beam-energy spread histogram is fitted by the function BES x with a binned log-likelihood fit with ROOT version 5.34.01 [21]. Figure 8a shows the best fit to the beamenergy spread with this model, and the resulting parameters are given in Table 1. The histogram contains 300 000 entries.
The width of the Gaussian σ and the boundaries of the beam-energy spread beta-distributions (x min , x max ) are fixed for all following fits. This assumes an existing knowledge of the beam-energy spread coming from the accelerator. Fixing these parameters can introduce a large systematic error, if they are not measured correctly.

Luminosity Weighted Beam-Energy Spread
The correlation between the particle energy and its position in the bunch causes a change in the effective beam-energy spread. The probability to radiate Beamstrahlung photons, and therefore the fraction of energy lost by particles, increases with the distance travelled in the electromagnetic field of the oncoming bunch. Figure 8b shows two energy distributions, one of the Peak-region, where both particles posses an energy of more than 99.5% of the nominal beam energy, and one of the Arms-region, where only one of the particle contains more than 99.5% of the nominal beam energy. Both histograms contain 300 000 entries.
The energy spread of the Peak-region is clearly flatter than the energy spread coming from the accelerator (Figure 8a). In the column Peak, Table 1 lists also the parameters ω a and ω b found by fitting Equation (14) to the distribution. For this fit, the limits and the Gaussian width are fixed. For the best fit both ω a and ω b are closer to zero, but still negative, i.e., there are still two maxima at the lower and upper end of the distribution. The peak at the lower end of the spectrum is reduced, because the particles in the tail are less likely to interact with particles that did not radiate Beamstrahlung.
The energy spread of the Arms-region, where one of the particles radiated Beamstrahlung, shows a large peak at the lower energy, and almost no peak at the upper end of the spectrum. This is also caused by the correlation between the energy and the position in the bunch. Particles in the tail are more likely to collide with a particle that already radiated Beamstrahlung, and therefore the peak at the lower edge of the beam-energy spread is enhanced. Likewise, only very few particles with the highest energy -located near the front of the bunch -interact with particles from the tail of the bunch, which leads to the disappearance of the peak at the highest energies. The beam-energy spread for the Armsregion is described by a beta-distribution for which ω b > 0 (see column Arms in Table 1).
The χ 2 /ndf becomes larger for the fits to the luminosity weighted beam-energy spreads than for the fit to the initial beam-energy spread. The chosen function cannot perfectly model the distributions, however, the fits are only used to check qualitatively if the model can represent the luminosity spectrum at this stage. In addition, as there was only a single input file available to run GUINEAPIG, the macro-particles Peak Energy spread in Fit + 99% C.L.
(b) Fig. 8 (a) The beam-energy spread distribution from the accelerator simulation [12] and the best fit to the beam-energy spread with Equation (14). (b) Energy spread after the simulation in GUINEAPIG. The distribution of the Peak requires that both particles have E > 0.995E Beam , and the distribution in the Arms requires that one particle has E < 0.995E Beam . The colour bands in both plots indicate the confidence interval at 99%.  (14) to the beam-energy spread from the accelerator simulation and to the beam-energy spread for two different regions of the luminosity spectrum.

Energy Spread
Peak Arms are re-used for luminosity events. This re-use means that the fluctuations in the number of entries are larger than what can be expected from the statistical uncertainties, which also increases the χ 2 /ndf value.

Beamstrahlung
Following Ohl's CIRCE model [22], the energy distribution of the particles after the emission of Beamstrahlung photons is modelled with a beta-distribution. Beamstrahlung will always reduce the energy of a particle, so that the random variate ∆ x Strahl would be between −1 and 0 (cf. Equation (10)). Beta-distributions are limited between 0 and 1, so that the function describing the Beamstrahlung effect is convoluted with the δ -distribution from Equation (9), which moves the range to 0 < x Strahl = 1 + ∆ x Strahl < 1 and no further variable transform is necessary for the probability density function.
The parameters of the beta-distributions used to describe the energy distribution due to Beamstrahlung are called η a and η b . The beta-distribution parameters must fulfil the conditions 0 < η a and −1 < η b < 0 for the distribution to fall towards x = 0 and rise towards x = 1.
Previous studies by Daniel Schulte have shown that the tail of the CLIC centre-of-mass energy distribution is better modelled by a sum of three beta-distributions. Therefore, the energy distributions from Beamstrahlung are initially fitted by linear combinations of N Beta incomplete beta-distributions and the constraint where p i are the respective fractions of the individual betadistribution contribution and [p] i = {η a i , η b i } the parameterset for each beta-distribution. The beta-distributions are limited with an upper limit of β Limit = 0.995. Above 0.955 the beam-energy spread is dominant and would have to be included for the fit. Figure 9 shows the fits with N Beta = 1, 2, 3 to the distribution of the particle energy. In Figure 9a the fit to the histogram is performed in the range of 0.0 < x < 0.995; It is visible that the function with three beta-distributions -with eight free parameters -shows a better agreement with the distribution than the other functions. As all the betadistributions cover the full range for the fit, there are large correlations between the parameters of different beta-distributions. Figure 9b shows the same fit of linear combinations with a range limited to 0.5 < x < 0.995; all three fit-functions overlap. Therefore, a single beta-distribution is enough to describe the particle energy between half and 99.5% of the beam-energy. For the Model, the Beamstrahlung is described by a single beta-distribution to reduce the number of free parameters. However, this will also limit the energy range in which our Model can be considered as valid.

The Model for the Full Luminosity Spectrum
The individual contributions discussed in the previous sections are now used to create the Model of the basic twodimensional luminosity spectrum. As was discussed in Section 3.1 the beam-energy spread BES x is described by a convolution of a Gaussian function g x and a beta-distribu- The beta-distribution for the beam-energy spread has a very narrow range. The particle energy distribution including the energy loss due to Beamstrahlung is described by a convolution of the beam-energy spread with an incomplete betadistribution with the upper limit of β Arm Limit = 0.9999, The upper limit is chosen to be close to 1, so that the convolution with beam-energy spread causes an overlap with the Peak-region (cf. Figure 10). To describe particle energy distributions only negligibly affected by the beam-energy spread, a beta-distribution with an upper limit of β Body Limit = 0.995 convoluted with a Gaussian function is used This upper limit separates the distribution from those more significantly affected by the beam-energy spread. This function is different from (17) due to the different ranges of the beta-distributions. As described in Equation (11), the distributions in the four different regions are described by the product of two functions, one for each particle. The explicit piecewise description shown in Equation (11), however, is replaced by the use of delta-distribution and implicit ranges of the individual functions. The Peak region is described by two pure beam-energy spread functions (Equations (14) or (17)) and delta-distributions to signify the absence of Beamstrahlung; the Arms are modelled by one beam-energy spread function and a delta distribution, and one function describing Beamstrahlung convoluted with the beam-energy spread (Equation (18)); the Body is described by two functions describing only the Beamstrahlung (Equation (19)).
with β Arm Limit = 0.9999, β Body Limit = 0.995. In addition, the constraint has to be fulfilled, which results in as required for a probability density function. The function given in Equation (20) will be used to describe the luminosity spectrum. The random variates according to the individual parts of Equation (20) are shown in Figure 10. Each summand of Equation (20) corresponds to one of the distributions. Due to the convolution of beam-energy spread and Beamstrahlung functions, the region around the nominal beam energies (x 1 ≈ 1, x 2 ≈ 1) is described by a superposition of individual contributions.

Reweighting Fit
The separate one-dimensional parts of the luminosity spectrum were fitted to the parts of the Model. Now the complete Model has to be fit to the two-dimensional spectrum. It is possible to fit Equation (20) to the basic luminosity spectrum. The convolutions with the δ -distribution can be performed explicitly. The other convolutions have to be performed numerically, because the convolution between the beta-distribution and the Gaussian function cannot be expressed in a closed form 4 .
For the implementation of the function the numerical convolutions are evaluated with the QAG 5 integration algorithm [23] interfaced via the GSLIntegrator from ROOT MathMore. The evaluation of the function takes about 160 seconds for the full range. A direct fit with the function, requiring multiple iterations, would be slow. The fitting procedure can be sped up by using a reweighting fit and by exploiting the fact that the random variates according to Equation (20) can also be described by the sum of the random variates of the individual functions [24].
The principle of the reweighting technique is shown in Figure 11. A χ 2 minimization, utilizing MINUIT [25] implemented in ROOT, is used to fit a sample of Model events to the GUINEAPIG sample. The procedure starts with the gen- 4 We have no formal proof of this statement. However, neither the integral of the Gaussian function (resulting in the error-function) nor the integral of the beta-distribution (yielding Gamma-functions) can be expressed in a closed form with a finite number of elementary functions. 5 is used to weight each event of the Monte Carlo distribution. For every set of parameter values [p] N a reweighted Monte Carlo distribution is obtained. The minimum χ 2 between the distribution of the Model and the distribution from GUINEA-PIG corresponds to the optimal parameter values matching the GUINEAPIG sample.
The χ 2 between the two histograms is calculated from the number of entries in bin j of the GUINEAPIG sample N j GP and its uncertainty σ j GP , the sum of the weights in bin j of the Model sample N j Model , and the uncertainty σ j Model , Minimizer: Weight every event i with its weight w i which are calculated from the event samples according to The χ 2 to be minimized is then calculated with where is a scaling factor that takes into account the difference in the sample sizes and the normalisation of the event weights due to the limited number of Model-events. The entire procedure has the advantage that only one sample of Modelevents is needed, contrary to traditional template-fit procedures that require generating new Monte Carlo samples for every parameter-set. However, by itself this requires even more evaluations of Equation (20) -one for every Monte Carlo event used in the generated sample -but the random variates according to this function can also be described by the sum of the random variates of the individual functions. The particle energy can be built up from the individual contributions where x Strahl is the random variate from the beta-distribution for the Beamstrahlung, x Spread the random variate from the beta-distribution for the beam-energy spread, and x G the random variate from the Gaussian-function of the beam-energy spread. Each random variate can be generated according to its probability density function. Equation (27) is connected to Equation (10): x Strahl = x E Beam + ∆ x Strahl , and x Spread + x G = ∆ x Spread . During the generation of events, the combination of functions is chosen according to the probability given by the parameters for each region p Peak/Arm1/Arm2/Body . The probability for a particle's energy in an event is given by the product of all individual probabilities and the product of the probabilities for the individual particles multiplied by the probability for the region gives the probability for the event Thus Equation (23) becomes and no numerical convolutions have to be calculated. The probability for obtaining energies x 1 and x 2 is not the same as the probability to obtain a specific group of variates, even if There are many combinations of the variates x Spread , x Strahl , and x G , which can lead to the same x 1 or x 2 . To estimate the probability for any given pair of energies L x 1 , x 2 the convolutions have to be performed either numerically or via Monte Carlo generation.

Application of the Reweighting Fit to Other Distributions
The reweighting fit is also used to fit the distributions after the inclusion of the Bhabha scattering, Initial and Final State Radiation, and detector resolutions. The individual events are passed through the Bhabha Monte Carlo generator and detector simulation, which can be understood as additional convolutions of the existing distribution. As can be seen in Equation (30), if the parameter governing one of the contributions does not change, the contribution does not affect the new weight. This enables the use of the reweighting fit also for the reconstruction of the spectrum from the Bhabha events, because the Bhabha scattering and detector resolutions D O k are not varied during the fit. A measured distribution f of Observables O k can be approximately written as where σ represents the centre-of-mass energy dependence of the Bhabha scattering and the observables, ISR E 1 , E 2 represents the probability for the energy distribution after Initial State Radiation, and FSR O 1 , O 2 , . . . represents the probability for the energy distribution after Final State Radiation. If the cross-section and detector resolutions are well enough known, the only difference between the measured and generated distributions is the luminosity spectrum. For this study the same Bhabha generator and detector simulations are used for both samples, so the additional effects are statistically the same. Any difference for the contributions can lead to a systematic error in the reconstruction of the luminosity spectrum

Equiprobability Binning
The χ 2 -fit requires binned histograms. To obtain an unbiased estimator of the compatibility in a χ 2 -fit, all bins should contain at least seven entries, and the number of events in all bins should be similar [26, p.304]. These requirements can be fulfilled when an equiprobability binning is generated based on the respective GUINEAPIG sample used in the fits. With equal-size bins either a large number of bins could be used -where most would contain very few or no entries and would have to be rejected for the χ 2 calculation -and the peak substructure could be resolved, or fewer bins with larger dimensions could be used, but then the peak could not be resolved. Therefore, the equiprobability binning can make better use of the available events.
Following the recipe of James [26, p.305], the events are first evenly separated along one axis, and then all events falling in the range on this first axis are again evenly separated in the second axis. If additional dimensions are used, the separation is repeated. In this way each bin has different dimensions along each axis, but the number of events per bin is constant.
For the fit to the basic luminosity spectrum, as discussed in Section 5.1, the distribution of the two particle energies is stored in a two-dimensional histogram. For the reconstruction of the spectrum from the Bhabha events the energy of the scattered electron and positron, and the relative centre-of-mass energy reconstructed from the acollinearity √ s acol/ √ s nom are filled into a three-dimensional equiprobability histograms. Figure 12 shows examples for a two-and threedimensional bin structure. It can be seen that around the nominal beam energies the size of the bins becomes smaller. Because the separation of events is done individually along each axis, the bin structures are not symmetric.

Uncertainty Estimation
In order to ensure that the Model is unbiased and consistent, a large number of Model vs. Model fits were performed with a varying number of bins. In each case, the procedure is as follows: two sets of events are created according to the Model of the basic luminosity spectrum. The samples are then used in the fit procedure described before. In each fit a cut on the centre-of-mass energy of √ s > 1.5 TeV is applied.
The pull distribution of every free parameter is obtained and fitted with a Gaussian function. The Model is unbiased, if the mean of every pull distribution is close to zero. The uncertainty is correctly estimated, if the pull width is compat- ible with unity. This is called the Normality condition [26, p.310]. For all parameters the pull distributions are independent of the binning. Most parameters are unbiased (i.e., the mean is zero). Exceptions are the η a Body1 and η a Body2 parameters, whose pulls are not normally distributed. These parameters describe the behaviour of the beta-distribution at the lower edge of the respective particle energy distribution. Therefore, the bias is caused by the cut on the centre-of-mass energy, which reduces the sensitivity to the lower energy part of the Body. The lower limit of these parameters is zero, which is often found by the minimizer instead of the nominal value. When the cut on the centre-of-mass energy is removed, the pulls are symmetric, and the parameters are correctly estimated. It is also found that the width of the Gaussian function is consistent with unity, so the uncertainties are correctly estimated by the minimization procedure.

Luminosity Spectrum Reconstruction
All ingredients for the reconstruction of the luminosity spectrum -the Model and the reweighting procedure -are now available.
The fits based on the basic luminosity spectrum (Section 5.1) are used to assess the similarity between the Model and the GUINEAPIG spectrum. The energies of the electron and positron pairs are filled into the two-dimensional equiprobability structure used in the reweighting fit. In the next step, the cross-section scaling, the Bhabha scattering, and the smearing for the detector-resolutions are applied and the reweighting fit is done with the observables defined in Section 2.4. The step-by-step inclusion of the intermediate steps and their impact on the reconstruction is detailed elsewhere [27].
The initial values of the parameters, used to generate the Model events, are given in Table 2. All the regions are chosen to start with a similar number of events (25%). The start-  ing ω parameters are taken from the fit to the beam-energy spread before the collisions ( Table 1). The other parameters are chosen arbitrarily in a way to cause a behaviour similar to the GUINEAPIG luminosity spectrum. The position of the two boundaries for the beam-energy spread (x min and x max ) are also taken from Table 1. Table 2 also lists the lower and upper bounds limiting the values for the minimizer.

Fit to the Basic Luminosity Spectrum
To verify that the Model can represent the basic luminosity spectrum from GUINEAPIG, the distribution of the initial particle energies are used in the χ 2 -fit. The data histogram is shown in Figure 2. The Monte Carlo sample is shown in Figure 10. The GUINEAPIG sample consists of 3 million events and the Model provides 10 million events. Fits are done with a binning varied from 50 × 50 bins to 300 × 300 bins in steps of 10 bins. Only events with √ s > 1.5 TeV are used in the fit. The cut is applied because the Model has limited validity range, and the events below half the nominal centre-of-mass energy would have a negative impact on the fit result.
As an example for the result of the reweighting fit, Figure 13 shows a small section of the histogram mapped onto one dimension and the pull distribution for all the bins before and after the fit. The data histogram has a constant number of events per bin, as designed by the equiprobability binning. The pull distribution after the fit converged is well centred around 0 with a width of 2.3. The width of the pull distribution is not equal to 1, because the χ 2 /ndf is larger than 1. This means that the Model is not completely identical to the GUINEAPIG distribution. Two of the differences are the limited number of beta-distributions used to model the tail of the spectrum (see Section 3.2), where deviations appear, and the differences in the peak of the spectrum (see Section 3.1), where a much larger number of parameters would be needed. As the χ 2 /ndf is not equal to unity, additional parameters would enable a better description of the spectrum.
For one fixed binning, the fit to the luminosity spectrum was done 198 times with the same GUINEAPIG sample and independent Model samples. All the parameter values vary within their uncertainties. Therefore, the Model and the fit procedure are consistent.

Fits to the Observables
The observables are defined in Section 2.4. Binnings 6 from 10 × 10 × 10 bins to 80 × 50 × 50 bins were used in the fits. The binning step is 5 bins for the relative centre-of-mass energy and 10 bins for the particle energies.
For the last step, the Bhabha events generated with the scaled luminosity spectrum are smeared with the detector resolutions as described in Section 2.4. The selection cuts had to be modified, and the cut is applied on the centreof-mass energy calculated from the smeared four-vectors √ s 4-vec > 1.5 TeV and in addition on the individual particle energies E 1 > 150 GeV and E 2 > 150 GeV. To recover Final State Radiation, the energy of all photons in a 3 • cone around an electron is summed up.

Discussion of the Results
For the two stages of the reconstruction multiple fits with different binnings were done. However, as the reconstructed spectra are fairly similar, only one reconstructed luminosity spectrum per stage is shown in detail. In addition, the parameter dependence on the number of bins is shown. For the fits to the basic and scaled luminosity spectrum the results with 100 × 100 bins are shown. For the reconstruction from the observables the fits with 40 × 50 × 50 bins are shown.
In Table 3 the χ 2 /ndf and parameters extracted by the selected fit stages are listed. The reconstructed parameters are far away from the initial values of the parameters (cf. Table 2), therefore the fit results are not artificially improved by using a good starting point. The final values of the beamenergy spread parameters ω are close to the values found by the one-dimensional fit to the beam-energy spread distributions detailed in Section 3.1. Because of the cut on the minimum centre-of-mass energy, the sensitivity on the lower Beamstrahlung parameter η a is lost, and both fits give a result of ≈ 0 with large uncertainty for these parameters. The reconstruction of the upper Beamstrahlung parameter η b is consistent.
The largest variation in the parameters is observed for the beam-spread parameters ω . This increase is mostly due to the detector effects. In total the uncertainty increases by a factor ten, and the values are significantly different.
There are significant correlations between the parameters. The largest correlations are between parameters from the same beta-distribution. The correlations also increase when the additional effects are taken into account. Some of the changes of the parameter values could, therefore, be due to increased correlations. Table 4 lists the fraction of events with a centre-of-mass energy larger than 0.99 √ s nom from GUINEAPIG and from selected fits of the different fit stages. The uncertainty of the GUINEAPIG value is the statistical uncertainty from one million events. The uncertainty for the fits is calculated from the uncertainty of the individual parameters and accounts for the correlation between them.
The difference of the fractions between GUINEAPIG and the Model is less than one percentage point. Given the size of the uncertainties the difference is significant. However, processes with lower cross-section will effectively use smaller  (24) and (26)) between the GUINEAPIG and Model samples before and after the re-weighting fit. A Gaussian function is fitted to the distribution of pulls after the fit. Table 3: The parameter values found in selected fits to the initial electron and positron energies (first rows) and to the observables (second rows). The details of the fits are given in the text. The basic luminosity spectrum from GUINEAPIG compared with the reconstructed basic luminosity spectra from the two fit stages for the selected fits are shown in Figures 14  and 16. For the ratios the green error bars show the statistical uncertainty for one million GUINEAPIG events and the barely visible red error bars show the uncertainty coming from the parameterisation.
In both cases the luminosity spectrum is reconstructed within 5% between 0.55 √ s nom and 0.995 √ s nom . Close to the peak, above 0.995 √ s nom , the beam-energy spread is the dominant effect and the difficulty of modelling this peak becomes visible. Still, this difference is seen only, when looking at small bin sizes (e.g., compare the bins around 1 in Figures 14d or 14e with Figure 14f). As Table 4 shows, the average fraction around the peak is reconstructed within 1 per-centage point. Improved parameterisations should be able to better describe and reconstruct the shape of the peak, at the cost of longer run-time for the fit.
Below 0.5 √ s nom , the Model is much more inconsistent with GUINEAPIG, but this is given by the design of this Model and the cut on the centre-of-mass energy applied for the fits.
Some of the reconstructed parameter values depend on the number of bins used in the fit. Figure 15 shows the dependence of the reconstructed parameters on the number of bins used in the fit. Fits with a binning of 50 × 50 bins to 300 × 300 bins with the same number of events were done. In Figure 15  with more bins the parameter p Peak rises, while the two parameters p Arm1 and p Arm2 fall, which is also visible in the correlation matrix and their correlation coefficient of about −0.4. Figure 17 shows the parameters obtained in the fit to the observables. The results are again sorted by decreasing χ 2 /ndf, which defines the Binning ID. In the figure the different markers give the number of bins used for the energy observables. As the χ 2 /ndf falls with increasing number of bins the larger the Binning ID the large is also the number of bins used for the relative centre-of-mass energy observable.
The parameter values depend much stronger on the number of bins. This is mostly due to the inclusion of the detector resolutions. Without a minimum number of bins the peak structure cannot be resolved, and the ω -parameters are completely different from the previous results and show large fluctuations in their values. If a large enough number of bins is used, the results are only a few sigma different from the previous fit results. The detector resolutions have a strong impact on resolving the structure of the luminosity peak.

Systematic Impact on Smuon Mass Measurement
There are significant differences between the reconstructed luminosity spectrum and the one from GUINEAPIG when looking at large event samples. Typical cross-sections for New Physics phenomena will be much smaller than that of Bhabha scattering, and the luminosity spectrum sampled for a specific process will therefore have larger statistical fluctuations, so that the difference between the reconstructed and actual spectrum might not be significant. To estimate the impact of the difference between GUINEAPIG and the reconstructed spectrum, the measurement of the smuon mass mμ± and neutralino mass mχ0 from smuon pair production is used. In this model the masses are mμ± = 1011 GeV and mχ0 = 340 GeV.
The smuon decays into a muon and a neutralino, so that the energy spectrum of the muons f (E µ ) can be used to extract the smuon and neutralino masses. The details of the analysis are described elsewhere [28], here only the parts directly concerning the systematic uncertainty from the luminosity spectrum are repeated. There are some differences in the treatment of the statistical uncertainty between the version of the fitting program used here, and the one used in the original paper.
In an ideal situation -with a single centre-of-mass energy √ s nom -the muon energy spectrum is a uniform distribution U E µ with the boundaries [29] The uniform distribution therefore depends on the smuon and neutralino masses. In reality, there is not a single centre-of-mass energy, and for every centre-of-mass energy the uniform distribution has different limits. Therefore, the measured muon-energy spectrum is affected by the basic luminosity spectrum, the Initial State Radiation, the cross-section, and the detector resolution D E µ . The luminosity spectrum L x , Initial State Radiation ISR x , and cross-section σμμ √ s can be combined into the number of events per centre-of-mass energy N eff (x). The Initial State Radiation and luminosity spectrum are convoluted and the resulting function is multiplied with the smuon-pair production cross-section and with the total integrated luminosity L int The Initial State Radiation function ISR describes the distribution of the energy after the radiation of initial state radiation. In this case, the distribution is obtained from the Monte Carlo generator used to generate the smuon events. It is here assumed to be independent of the nominal centre-ofmass energy. The function to fit the muon energy spectrum is then the convolution of the uniform energy spectrum with the detector resolution weighted by the respective number of events (34) Figure 18a shows the background-subtracted signal sample and an example fit with Equation (34). To estimate the impact of the reconstruction, the fit results when the luminosity spectrum is taken directly from GUINEAPIG are compared with those, when the spectrum is coming from the reconstruction. The masses extracted from the fit with Equation (34) become a function of the parameters p from the spectrum reconstruction m = m(p) with the luminosity spectrum reconstructed from the Bhabha events. To estimate the systematic uncertainty due to the reconstruction of the spectrum, the fit is performed with the nominal set of parameters p and with each parameter p i increased or decreased by half of a standard deviation σ p i The systematic uncertainty on the fitted value is then given by i , and the correlation matrix C. Table 5 lists the smuon and neutralino masses from the fit when the luminosity spectrum in equation (33) is directly taken from GUINEAPIG and when the luminosity spectrum is obtained from the reconstruction with the observables with the scaled luminosity spectrum and detector resolutions with a binning of 50 × 40 × 40 bins. The difference in the reconstructed masses for these two luminosity spectra is smaller than the statistical uncertainty. However, as the reconstructed luminosity spectrum shows a dependence on the binning, so do the reconstructed masses. Figure 18b shows the reconstructed masses for the spectra reconstructed with different binnings. There is a dependence of the reconstructed masses on the number of bins, but the spread of the reconstructed masses is smaller than the statistical uncertainty (cf. Table 5).
As the difference between the obtained masses and the spread of masses is smaller than the statistical uncertainty, the reconstruction of the luminosity spectrum does not introduce a significant bias compared with the statistical uncertainty. The systematic uncertainty due to the luminosity spectrum reconstruction is also much smaller than the statistical uncertainty, so that the total uncertainty on the reconstructed mass is not increased significantly.

Summary, Conclusions, and Outlook
A framework has been developed for the reconstruction of the basic luminosity spectrum at future linear colliders. The spectrum can be reconstructed from Bhabha events measured with the tracking detectors and calorimeters. All important effects were included: the luminosity spectrum from beam-beam simulations -including the non-Gaussian CLIC beam-energy spread -the √ s -dependence of the Bhabha cross-section, Initial and Final State Radiation, and the detector resolutions.
The Model of the 3 TeV CLIC luminosity spectrum, required for the reweighting fit, has some limitations. For technical reasons the energy range to describe the tail of the Beamstrahlung is limited to √ s > 1500 GeV, and the peculiar beam-energy spread cannot be modelled precisely with few parameters. The reweighting fit itself does not impair the reconstructed spectrum. The differences between GUIN-EAPIG and the reconstructed spectrum do not significantly change between the fit to the basic luminosity spectrum and the fit to the observables with the scaled luminosity spectrum and including detector resolutions. With an improved model, and increased processing power, an improved reconstruction of the CLIC 3 TeV spectrum should be possible.
The fraction of events above 99% of the nominal centreof-mass energy is reconstructed within 1 percentage point. The centre-of-mass energy distribution is reconstructed to better than 5% between the nominal and about half the nominal centre-of-mass energy, the validity limit of our Model. These results are obtained regardless of the included level of details, so that one can conclude that the limitations of the Model cause most of the discrepancies to the simulated spectrum, and if a better model is used, the discrepancies should be reduced.
To estimate the systematic impact on other physics measurements, the reconstructed spectrum was used in the study of smuon decays, one of the CLIC 3 TeV benchmark processes. The reconstructed spectrum does not induce a significant bias on the measured mass, nor does it cause a significant systematic uncertainty. The systematic uncertainty from the spectrum reconstruction is two orders of magnitude smaller than the statistical uncertainty.
The spectrum is well enough reconstructed for the chosen physics channel. In this case a good reconstruction of the tail of the spectrum is tested. The reconstruction of the peak is less important, because the process is far above threshold and the cross-section does not change significantly over the peak region. More work is needed to evaluate and possibly improve the reconstruction of the peak.

Outlook
The framework can also be applied for the reconstruction of the luminosity spectrum at other centre-of-mass energies and linear electron-positron colliders than CLIC. Depending on the beam-energy spread and the demanded range of the reconstruction, the Model has to be adapted, but this will not increase the computational complexity of the reconstruction. The energy range of the current Model can be increased by replacing the single Beamstrahlung beta-distributions by linear combination of beta-distributions. Improving the description of the beam-energy spread is less obvious without a large increase in the number of parameters.
The boundaries of the beam-energy spread -the parameters x min and x max -were fixed during the fit. It should be evaluated how much the measurement is affected, when these parameter values differ from those of the beam-energy spread. It should also be tried to vary the boundaries of the beam-energy spread during the re-weighting fit. For varying these parameters during the reweighting fit the initial samples have to be produced with overlapping regions. For example, the peak region would be produced with an x min smaller than the upper limit of the arm or body regions. During the re-weighting the value for x min or x max would be given by the minimizer, and events in the peak below x min or above x max would be dropped, as would events in the arm or body above their respective upper limit.
The observables from the Bhabha events can also be exchanged for other suitable choices, always keeping the detector resolutions in mind. The impact of the detector resolutions on the reconstructed spectrum can be easily studied by changing the resolutions used in the four-vector smearing. The same detector resolutions and Bhabha generator were used for the GUINEAPIG and Model events. Differences in the predicted detector resolution and Bhabha scattering to the actual events can introduce systematic errors into the reconstruction. These effects could be studied by varying the detector resolutions or the Bhabha cross-section independently for the two samples used in the fit.
Only Bhabha events -and no other physics processes -were considered. It should be checked if multi-peripheral two-photon events, in which the spectator electrons scatter at large angles, are a background.
As the luminosity spectrum depends on the accelerator, the impact of possible variations of the beam parameter on the reconstruction of the luminosity spectrum should be studied with realistic variations of the beam parameters.