Vincent Favaudon
email: vincent.favaudon@curie.fr

Jean-Michel Lentz
email: lentz.jeanmichel@gmail.com

Sophie Heinrich
email: sophie.heinrich@curie.fr

Annalisa Patriarca
email: annalisa.patriarca@curie.fr

Ludovic De Marzi
email: ludovic.demarzi@curie.fr

Charles Fouillade

Marie Dutreix
email: marie.dutreix@curie.fr

Ludovic De Marzi

Time-resolved dosimetry of pulsed electron beams in very high dose-rate, flash irradiation for radiotherapy preclinical studies

Keywords:

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Five years ago Favaudon et al. [START_REF] Favaudon | Ultrahigh doserate FLASH irradiation increases the differential response between normal and tumor tissue in mice[END_REF] reported for the first time that short bursts of relativistic electrons at very high dose-rate allowed sparing C57BL/6J mice from radio-induced lung fibrosis at doses that inevitably trigger the development of the disease in mice exposed to conventional dose-rate irradiation (CONV). This methodology, named FLASH, also spared neural stem cells and spatial memory following total brain irradiation [START_REF] Montay-Gruel | Irradiation in a FLASH: unique sparing of memory in mice after whole brain irradiation with dose rates above 100Gy/s[END_REF]. In contrast FLASH was as efficient as CONV with tumor control as an endpoint [START_REF] Favaudon | Ultrahigh doserate FLASH irradiation increases the differential response between normal and tumor tissue in mice[END_REF]. Therefore, FLASH dramatically increases the differential outcome between normal tissues and tumors at the advantage of the former. Protection against normal tissue injury without loss of the tumor efficiency has been confirmed in large animals, thus evoking a promising potential of the FLASH methodology in anticancer radiotherapy [START_REF] Vozenin | The advantage of Flash radiotherapy confirmed in mini-pig and catcancer patients[END_REF]. FLASH might actually widen the therapeutic window of radiotherapy in hypofractionated protocols involving large doses per fraction to enable dose escalation or reduction of late complications. It might also improve targeting mobile tumors, for instance in respiration-gated radiotherapy in patients treated for lung cancer.

Were FLASH to be used in the clinic routine, it would involve major changes in the design and control of the particle accelerators and the dosimetric techniques medical physicists are accustomed to, irrespective of whether low megavoltage electrons, very high energy electrons (VHEE), photons or even protons would be used. In particular, due to ion recombination, space charge and ion drift effects, ionisation chambers may not provide a proportional response with high radiation doses per pulse even if strategies were recently attempted to correct for it [START_REF] Petersson | High dose-per-pulse electron beam dosimetry -A model to correct for the ion recombination in the Advanced Markus ionization chamber[END_REF]. Indeed, any FLASH session can be achieved using either a single pulse with a high peak current, typically 200 mA allowing to deliver as much as 40 Gy in a 1 µs pulse, or a limited number of submicrosecond pulses delivered in sequence at 5-10 ms interval in such a way that the total beam-on duration should remain below 100 ms [START_REF] Montay-Gruel | Irradiation in a FLASH: unique sparing of memory in mice after whole brain irradiation with dose rates above 100Gy/s[END_REF].

The sensors and oscilloscopes used for beam monitoring and dosimetry under such conditions must offer a large bandwith and enough memory to provide collection, storage and display of the complete pulse sequences. In the following we provide a description of the linear electron accelerator (Kinetron) used for FLASH studies, and give an account of the available techniques that allow characterization of the beam profile and properties to provide a time-resolved, proportional measurement of dose delivery. Application to the calibration of non-linear systems such as chromogenic films for postponed dose control, is also described.

Numerical Monte-Carlo simulations to model the relative dose distribution in mouse organs irradiated using the Kinetron at both conventional and very high dose-rate, have been published elsewhere [START_REF] Lansonneur | Simulation and experimental validation of a prototype electron beam linear accelerator for preclinical studies[END_REF].

MATERIALS AND METHODS

-Accelerator design and operation

2.1.1-General features

The S-band linear electron accelerator (LINAC) that was used for the original FLASH studies [START_REF] Favaudon | Ultrahigh doserate FLASH irradiation increases the differential response between normal and tumor tissue in mice[END_REF], named Kinetron was designed and built by CGR-MeV Co. (Buc, France) in 1987 with a view to investigate the kinetics of electron transfer and radical reactions in macromolecules at the submicrosecond time-scale [START_REF] Favaudon | CO2-radical induced cleavage of disulfide bonds in proteins. A gamma-ray and pulse radiolysis mechanistic investigation[END_REF]. The system (Figure 1) was adapted to preclinical studies by EuroMeV (Buc, France) and PMB-ALCEN (Peynier, France). It incorporates an optical stage synchronized with electron pulses allowing time-resolved chemical dosimetry in water [START_REF] Favaudon | CO2-radical induced cleavage of disulfide bonds in proteins. A gamma-ray and pulse radiolysis mechanistic investigation[END_REF]. The Kinetron accomodates optical cuvettes, cells in culture flasks or animals including cat patients [START_REF] Vozenin | The advantage of Flash radiotherapy confirmed in mini-pig and catcancer patients[END_REF] and can deliver a wide range of pulse sequences, from a single submicrosecond pulse to repeated pulses at a known frequency, or split-dose irradiation with an interval between pulses ranging from 5 ms to a few minutes [START_REF] Ponette | Hyperfast, early cell response to ionizing radiation[END_REF].

The Kinetron is powered by a 5.0 MW thyratron/magnetron stage generating RF pulses of fixed duration (3.2 µs) at an adjustable repeat frequency, and is fitted with a thermoionic triode electron gun. Precise adjustment (1 V resolution) of the voltage applied to the triode grid allows a full control of the emitted current and pulse width in the FLASH mode, i. e., at a high beam current. In addition the Kinetron incorporates an automatic tuning of the magnetron RF frequency (AFC), and an adjustment of the kilovoltage applied to the magnetron (DQing).

The Kinetron is operated by a computer through an interface programmed with LabView™ (National Instruments, Austin, Texas), allowing a full control of the DQing, triggering, delays, pulse width and sequences, peak current and repeat frequency (Table 1). The dose per pulse at the target was set by adjusting the peak current, the pulse width and the distance to the electron output window. The pulse to pulse repeatibility at 100 mA peak current, was better than ± 1.2% in terms of delivered dose after the LINAC had reached a thermal equilibrium.

The Kinetron was also used to perform irradiation at dose-rates in the same range as conventional radiotherapy, i. e., at a few Gy per minute. In that case the peak current was set to 80 µA to provide a dose per pulse of about 1 mGy at 1 m from the source.

-Monitoring electron pulses at high beam current

The set-up for irradiation is shown in Figure 2 The software operating the Control-Command computer and the NI PXI-1033 digitizer was originally designed by EuroMeV, and allows real-time control of the LINAC power stages as well as automated collection, storage and integration of the charges collected by the primary collimator.

-Optical system

The Kinetron facilities involves an optical system (Figure 3) synchronized with electron pulses allowing time-resolved analysis of transmittance changes in liquid samples [START_REF] Favaudon | CO2-radical induced cleavage of disulfide bonds in proteins. A gamma-ray and pulse radiolysis mechanistic investigation[END_REF].

Signals from the photomultiplier tube were stored on the Nicolet oscilloscope. The risetime of the system was about 20 ns with an RC coupling of 50 Ω.The cuvettes containing the solutions were maintained in an holder positioned in the electron beam, in such a way that the electron and light beams penetrated the cuvette at 90° from each other in the horizontal plane. The cuvettes had 1.25 mm thick optically flat windows, with path lengths of 2 and 10 mm along the electron and light beams, respectively. The water-equivalent depth at the center of the irradiated liquid layer, taking into account the density of silica (d = 2.2), was 3.75 mm. Note that the electronic equilibrium is not reached at this depth and that a correction has to be made as long as calibration with respect to the Markus ionisation chamber or irradiation of thick samples, is addressed.

-Time-resolved dosimetry in water based on the methyl viologen radical

300 µl aliquots of a 3 mM solution of methyl viologen (1,1′-dimethyl-4,4′-bipyridinium dichloride) in 100 mM sodium formate (Normatom grade, VWR Chemicals, Fontenay-sous-Bois, France) made in ultra-pure water (Direct-Q3, Merck-Millipore, Molsheim, France) were introduced into fused silica cuvettes (see § 2.1.3), each fitted with a sealed tubing and gastight septum inlet adapter as described [START_REF] Favaudon | CO2-radical induced cleavage of disulfide bonds in proteins. A gamma-ray and pulse radiolysis mechanistic investigation[END_REF]. Pure (99.99%) N2O (L'Air Liquide-Alphagaz, Mitry-Mory, France) was saturated with water and introduced through Ø 0.8 mm stainless steel tubings directly into the liquid layer, allowing gentle bubbling (1 h in dim light). Any contamination by oxygen was carefully avoided at all steps. Under these conditions, the OH • , H • and e radicals generated from the onset of water radiolysis are converted into the methyl viologen radical cation MV •+ :

a q - H2O  Radiation → OH • + H • + e -aq
Water radiolysis e -aq + N2O + H2O → OH • + OH -+ N2

Conversion of e -aq into OH •

H • + N2O → OH • + N2 Conversion of H • into OH • OH • + HCOO -→ • COO -+ H2O Conversion of OH • into • COO - • COO-+ MV 2+ → CO2 + MV •+ Generation of the MV •-radical
Due to the large difference in the redox potentials of the HCOO -/ • COO -(E°' = -1.85 V) and MV 2+ /MV •+ couples (E° = -0.70 V) the radiolytic yield for the formation of the MV •+ radical is the same as that for • COO -, namely, G = 0.625 µmol.J -1 . MV •+ presents two aborption bands at 393 and 603 nm [START_REF] Das | Reactions of methyl viologen dication (MV2+) with H atoms in aqueous solution: mechanism derived from pulse radiolysis measurements and ab initio MO calculations[END_REF]. We used the 603 nm band (ε603 = 13.3 mM -1 .cm -1 ) for dosimetric studies.

-Gafchromic™ EBT3 films dosimetry

The strategy used for the intercomparison of the response of Gafchromic™ EBT3 films (Ashland Advanced Materials, Bridgewater, New Jersey) to the FLASH vs. CONV doserates, involved 7 steps.

To calibrate EBT3 films at CONV dose-rate we compared the response of a PTW (Freiburg, Germany) Advanced Markus ionisation chamber type 34045, to that of the methyl viologen dosimeter. The accelerator was set at 80 µA beam current and the following experiments were performed. (i) The Markus chamber was mounted in the PMMA stack (d = 1.17) provided by PTW, placed at the center of and perpendicularly to the electron beam at 798 mm from the source and exposed to 1,000 pulses at a repeat frequency of 100 Hz.

Layers of PMMA were removed until leaving the chamber bare in air within the PMMA slab.

The electronic equilibrium was reached at 7 mm PMMA (Figure 4A). Note that the PMMAequivalent depth relative to the center of the methyl viologen cuvette is 3.2 mm. (ii) The Markus chamber was set at electronic equilibrium in the solid water phantom at 500 mm from the electron output window. The dose per pulse at this distance was about 7 mGy. No change in the response of the Markus chamber was observed as the repeat frequency was increased from 10 to 250 Hz. Thus, a 4 ms delay between pulses is sufficient to collect to the electrodes 100% of the charges generated by 7 mGy pulses, consistent with the ion collection time (22 µs) guaranted by the manufacturer. (iii) The Markus chamber was moved along the axis of the electron beam, from 800 to 220 mm from the source, in such a way that the dose-rate varied between 7x10 2 and 10 4 Gy.s -1 during 3.2 µs pulses. There was no evidence of saturation at the maximum dose-rate used. (iv) The doses received at the Markus chamber were compared to those determined from the methyl viologen dosimeter at the same position, taking into account in-depth differences as in § 3.1.3. The formation of the MV •+ radical was measured by online optical detection at 603 nm, with correction to take into account the instability of the MV •+ radical as the length of radiation exposure was > 10 s (see § 3.2). (v) Pieces of EBT3 films were irradiated at electronic equilibrium in contact with the Markus chamber and exposed to dose escalation, from 1 to 40 Gy. The optical density in the films was measured by spectrophotometry (Figure 4B) and data blinded to avoid any tester's effect.

For calibration of EBT3 films at very high dose-rate, the beam current was set at 200 mA to provide 2x10 7 Gy.s -1 during the pulse at 220 mm from the electron output window. DQing was adjusted so as to obtain the same in-depth profile as at CONV dose-rate. The following experiments were performed by varying the pulse width, the number of pulses (repeat frequency 100 Hz) or the distance to the electron output window. (vi) The formation of the MV •+ radical was followed by online optical detection and compared to the signal integral measured at the primary collimator. The beam-on time was below 500 ms, in such a way that no correction for the secondary decay of the MV •+ radical (see § 3.2) was necessary. (vii) The experiment was repeated by substituting EBT3 films for the methyl viologen cuvette using a PMMA stack providing the same water-equivalent depth (3.75 mm). Data were blinded.

24-h after irradiation the films were digitized using either an Epson Expression 11000XL, 48-bit color depth flatbed scanner (Suwa, Japan) and FilmQA Pro software following the OneScan protocol (Ashland Advanced Materials), and/or by recording the optical absorption spectrum of the films using a Hitachi (Tokyo, Japan) U-2000 spectrophotometer.

-Lynx ™ scintillation detector -Application to the monitoring of beam profiles.

The Lynx ™ detector (FIMEL, Fontenay-aux-Roses, France) is a planar, 2-D system for charged particles dosimetry based on a high resolution scintillating screen coupled with a charge-coupled device (CCD) camera. Indeed, plastic scintillators are energy independent at the energies found in clinical electron beams and photon beams and poorly sensitive to the angle of incidence of radiation, making them good candidates for relative dosimetry in high energy charged particle beams [START_REF] Beaulieu | Review of plastic and liquid scintillation dosimetry for photon, electron, and proton therapy[END_REF]. As an example, full dosimetric characterization of the Lynx ™ detector has been described for scanned ion beams [START_REF] Russo | Characterization of a commercial scintillation detector for 2-D dosimetry in scanned proton and carbon ion beams[END_REF].

The detector has an active surface area of 300 x 300 mm with an inherent spatial resolution of 0.5 mm. The entrance face of the detector is composed of 2.4 mm glass fiber and 70 µm copper layers. The scintillating screen (0.5 mm thickness gadolinium-based plastic material) converts the energy lost by the radiation into photons which are reflected by a mirror and collected by the photodetectors of a CCD camera (a 540 nm wavelength filtration is obtained limiting the spectral analysis). A variable aperture collimator system (iris) allows modifying the amount of light reaching the CCD camera. The detector signal from the CCD camera is digitalized with 10-bit depth. Dedicated software Lynx 2D ™ is used for image acquisition, raw data correction for image geometrical distortion and background noise.

Single frame images with adjusted exposure time can be collected to obtain 2-D dose distributions from the scintillation screen.

-Small-field dosimetry using Mallet-Cerenkov light emission

To investigate the feasibility and linearity of radiation dosimetry in small fields using Cerenkov light emission, we designed a Ø 14x80 mm duraluminum tube in which were assembled, at one end an optical mount comprising a Ø 10x14 mm fused silica cylinder (n = 1.46, d = 2.2) with optically polished faces, a biconvex adaptator lens (f = 10 mm) also made of fused silica and held by washers, and an aluminum foil maintained by a drilled cap to protect from ambient light; at the other end a bundle of seven 1 mm PMMA optical fibers, 12m long inserted into a duraluminum hollow cylinder whose external diameter was adjusted to 10 be moved along the master tube to recover the maximum amount of Cerenkov light at the optimal distance from the focal plane of the lens (Figure 5). The length of the fused silica cylinder was calculated for all incident electrons to be stopped within the silica core. The light emerging from the terminal end of the optical fiber bundle was received at a Hamamatsu Photonics (Massy, France) R928 photomultiplier optimized for the detection of nanosecond transients. The probe was positioned along the axis of the electron beam, at a known distance from the source and exposed to single or sequential pulses of varying current, width and energy.

-Data analysis

Data were analyzed using Kaleidagraph software (Synergy software, Reading, Pennsylvania)

to fit time-dependent processes to linear or non-linear mathematical equations according to the Gauss-Levenberg-Marquardt algorithm [START_REF] Press | Numerical Recipes[END_REF]; unless otherwise stated, the results were expressed as mean ± standard deviation. The goodness of fit relative to the functions y = ax or y = ax + b, was appreciated using a Spearman's test with StatEL software (AD Software, Paris).

RESULTS

-Beam monitoring

3.1.1-Fluence measurement at high beam current

Screenshots and online computer analysis of the signal generated by electron pulses in the primary graphite tore (Figure 2, item #4) are shown in Figure 6. These experiments establish the proportionality of the charge collected at the tore with the pulse width as well as the number of pulses fired in sequence.

-Depth-dose distribution

To mimick the dose distribution in the depth of irradiated targets, stacks of Gafchromic™ EBT3 films (see § 2.3 and 3.3) tightly packed between 3.0 mm thick lucite slabs (d = 1.36)

were set at 500 mm from the source and exposed to 10 Gy (dose in air) at various DQing voltage. The doses were determined from the optical absorption at 583 nm in the films as described in § 3.3. Figure 7 shows the build-up and attenuation of the dose as a function of the depth in the stack, and how the energy of electrons can be adjusted by modifying the DQing to compensate for the loss of energy at high current. This control was performed routinely along in vivo studies [START_REF] Favaudon | Ultrahigh doserate FLASH irradiation increases the differential response between normal and tumor tissue in mice[END_REF] to make sure that the depth dose distribution in tissues, hence the incident energy was the same in both radiation modalities.

-Issues due to small aperture collimators

In preclinical assays involving small laboratory animals, the use of a secondary collimator to limit the radiation field to the organs of interest, is a necessity. Volume averaging effects, or the lack of secondary electron equilibrium in the vicinity of such collimators may affect the penumbra and the electron spectrum or the divergence of the beam in a complex way and requires a precise determination of the dose distribution in each individual mounting.

A 22-mm thick graphite plate, sufficient to stop 5 MeV electrons was drilled with a rectangular 18x26 mm hole (HxL) in its center, and set perpendicularly to the beam axis. The plate could be moved back and forth by translation along the beam axis. 40x40 mm cuts of Gafchromic™ EBT3 film were put with smooth adhesive tape over the surface of a PMMA block positioned along the beam axis at 500 mm (fixed position) from the electron exit window of the Kinetron. The Kinetron was set to deliver 4 pulses cumulating 9.6 Gy at this position. EBT3 films were revealed 24 h post-irradiation by scanner or spectrophotometry.

The dose received at the films decreased by as much as 31% as the distance to the rear face of the collimator increased whilst the lateral penumbra increased with the distance (Figure 8). This air gap factor is attributed to the loss of side-scatter equilibrium, collimator scatter and other scatter effects. Beyond 35 mm the dose received at the films did not change any more and did not change either as the collimator was removed. In other words, beyond 35 mm from the collimator the dose was the same as that received from the direct beam in the absence of secondary collimator. For this reason the anterior face of irradiated targets should always be set at a suitable distance from the secondary collimator, in that case 35 mm, depending on the energy of electrons and the size of the collimator aperture.

-Time-resolved dosimetry in aqueous solution based on methyl viologen

The dosimetry of electron pulses was carried out by time-resolved, optical detection at 603 nm of the blue methyl viologen radical cation MV •+ in N2O-saturated formate buffer (see §

2.2).

The formation of the MV •+ radical was complete within the rise time of the photomultiplieroscilloscope stage (Figure 9A) and was stable for enough time to allow for determination of the dose delivered during successive pulses fired in sequence (Figure 9B). We also used the MV •+ radical to determine the inverse-square law fitting the intensity of the radiation as a function of the distance to the electron output window (Figure 9C).

However, following its formation in the submicrosecond time range the MV •+ radical decayed slowly with loss of the absorption at 603 nm in the course of a biphasic process (Figure 10A) that fitted a bi-exponential equation of the form:

(1) A(t) = a1.exp(-.t) + a2.exp(-k2.t)

where A(t) is the absorbance at time t, and a1 and a2 the amount of the 603 nm absorbance involved in the first and second step of the reaction with rate constants k1 and k2, respectively. Non-linear regression drawn for best fit to experimental data gave a1 = 0.0673 ± 0.0013, k1 = 1.071 ± 0.049 x 10 -2 s -1 (0.642 ± 0.029 min -1 ); a2 = 0.193 ± 0.001, k2 = 4.478 ± 0.582 x 10 -5 s -1 (0.00269 ± 0.00035 min -1 ); a1/(a1+a2) = 0.259, a2/(a1+a2) = 0.762. The decay was half-over in 8802 ± 78 s (146.7 ± 1.3 min). Note that (a1+a2) = D * G * ε where D is the dose, G the radiolytic yield for the formation of the MV •+ radical and ε its molar extinction coefficient at the wavelength of interest (see § 2.2).

Let us consider a situation in which the cuvette containing the methyl viologen solution is exposed to N pulses of elementary dose d each, delivered in sequence at a repeat frequency ν. The interval between two consecutive pulses is 1/ν, the time t at which n pulses have been fired is t = n/ν, the time at the end of the pulse sequence is T = N/ν and the time elapsed between any pulse and the last pulse is (T-t) = (N-n)/ν. By substituting these variables into equation ( 1) and summing over the number of pulses n it comes:

(2) 2) has no analytical solution but can be solved through a parametric method for any given value of N and ν. An example is given in Figure 10B.

A(N) = d * G * ε * { a1 (a1+ a2) exp[-k1 n=1 n=N ∑ * (N-n) ν ] + a2 (a1+ a2) exp[-k2 n=1 n=N ∑ * (N-n) ν ] } Equation (
In conclusion, the methyl viologen dosimeter allows time-resolved determination of the dose received at any position in the electron beam at high electron fluence in temporal coincidence with the signal received from the primary collimator. However, owing to the evolution of the MV •+ radical over extended periods of time, the use of the methyl viologen dosimeter is not recommended as long as the irradiation protocol consists of a protracted exposure to radiation, unless the decay of the MV •+ radical is taken into account according to eq. ( 2).

-Calibration of Gafchromic ™ EBT3 films with regard to the dose-rate

For preclinical investigations Jaccard et al. [START_REF] Jaccard | High dose-per-pulse electron beam dosimetry: usability and dose-rate independence of EBT3 Gafchromic films[END_REF] used Gafchromic™ EBT3 films set in contact with the animals and calibrated at low dose-rate with reference to a PTW (Freiburg, Germany) Markus Advanced ionisation chamber. However, calibration of EBT3 films at very high dose-rate was not done for lack of a system allowing time-resolved measurement of the incident dose.

We took advantage of the optical facility described in Figure 3 to calibrate Gafchromic™ EBT3 films with reference to the methyl viologen dosimeter following the strategy described in § 2.3. The optical density of the films exposed to any radiation dose in the CONV and FLASH mode were indistinguishable from each other within ± 1.84% precision (Figure 11).

Within these limits the photometric response of Gafchromic™ EBT3 films is independent on the dose-rate for up to 2x10 7 Gy.s -1 during the pulse. This experiment was repeated 4 times over a period of 3 years with different batches of EBT3 films. The sensitivity of the films to radiation or aging varied between batches, but the response of each batch was independent on the dose rate and the dosimetric precision was in the range ± 2% throughout.

-Monitoring beam profiles using the Lynx ™ scintillation detector

The transverse dose distribution of the electron beam, as well as its characteristics in terms of linearity with the pulse width and dependence on the grid voltage were measured with the Lynx ™ mounted at various positions in front of the emission window, perpendicularly to the beam axis. The "in-pulse" (peak) dose-rates were set at 0.4-3.5 10 6 Gy.s -1 throughout. Data analysis was performed in MATLAB ™ Version R2012a.

The results (Figure 12) show that the response in dose was linear in the range explored.

The FWHM value measured at 600 mm from the electron output window was 135 mm. The symmetry of the beam was assessed as the FWHM variation at each angle on a transverse beam profile and was within 1% for all tested conditions (high and low dose-rate).

The Lynx ™ system demonstrated an excellent linearity of the dose-response relationship with good reproducibility and stability. However, the CCD camera of the Lynx ™ suffered from a relatively narrow dynamic range. Substituting the CCD camera for a multichannel photomultiplier, and miniaturisation might allow one to overcome these pitfalls and allow precise small-field dosimetry.

-Small-field dosimetry using Mallet-Cerenkov light emission

It has long been known from the pioneering works of Mallet [START_REF] Mallet | Étude spectrale de la luminescence de l'eau et du sulfure de carbone soumis au rayonnement gamma[END_REF] and Cerenkov [START_REF] Cerenkov | The visible glow of pure liquid under the effect of hard gamma-rays[END_REF][START_REF] Jelley | Cerenkov radiation and its applications[END_REF][START_REF] Glaser | Applications of Cherenkov light emission for dosimetry in radiation therapy[END_REF] that a cone of UV light is emitted along the path of relativistic charged particles as those travel faster than the local speed of light in a dielectric medium. The angle θ at which the radiation is emitted relative to the particle's direction is given by the equation:

(3) cos θ = 1/βn where θ is the opening angle of the cone, n the refractive index of the medium and

β = v/c
where v is the velocity of the particle and c the speed of light in vacuum. The intensity of the Cerenkov radiation produced by a given electron fluence is proportional to (1 -1/β 2 n 2 ).

Accordingly the efficiency for Cerenkov light emission is nearly independent of the electron energy above 2 MeV but is expected to decay abruptly below 1 MeV [START_REF] Beddar | Cerenkov light generated in optical fibers and other light pipes irradiated by electron-beams[END_REF]. On the other hand the spectral distribution of Cerenkov radiation is broad, covering the entire visible range, with a maximum in the blue region.

In a first experiment the Cerenkov probe (Figure 5) was set at 850 mm from the electron output window of the Kinetron and exposed to 1.00 µs pulses, 0.47 Gy each at two values of the DQing (Table 1) such that the nominal energy of electrons was either 3.9 or 5.0 MeV.

Oscillocope recordings of the signal received at the photomultiplier through a 50 Ω load resistor showed a square signal with rapid oscillations of 35 ns period at the onset and the end of the pulse (Figure 13A). The period of these oscillations was equal to that expected for the light generated at the silica probe to transit through the optical fiber to the photomultiplier tube, taking into account the refraction index of the material constituting the fiber, suggesting that oscillations originate from multiple reflexion of Cerenkov light at the air-fiber interfaces.

Due to decoherence, the phenomenon faded away in about 500 ns. It did not affect the integral of the signal under the curve (∑V * t) and could be easily suppressed by applying a RC filter or an analog low-pass active filter. This notwithstanding, the integral was strictly proportional to the dose received at the silica probe. Interestingly, the integral was lower as the mean energy of electrons was decreased from 5.0 to 3.9 MeV as expected from the fact that reduced penetration of electrons in the silica block produces a decrease of the amount of light focused by the collecting lens over the fiber bundle. Thus, with a proper mounting the Cerenkov light could be used to detect occasional energy drops of the LINAC RF power source.

In a second experiment the probe was set at 300 mm from the electron output window.

The beam current was adjusted to provide 4.9 Gy in a 1 µs pulse at the silica block. The load resistor at the photomultiplier was set at 50 Ω to provide fast response. Single pulses were fired with stepwise increase of the pulse width from 0.10 to 2.20 µs. The beam current measured at the primary collimator and the Cerenkov signal measured at the photomultiplier, were recorded in parallel. The result (Figure 13B) shows that the integral of the Cerenkov signal increases with the pulse width, in exactly the same way as the charge collected at the collimator, and with no evidence of saturation.

In a third experiment the beam current was lowered down to 100 µA. The probe was set at either 300 or 806 mm to receive either 12.8 or 1.77 mGy per pulse, respectively. The gain of the system was sufficient to detect the beam current and pulse profile with a load resistor set at 50 Ω (Figure13C). The signal integral over the pulse increased proportionally to the ohmic load applied to the photomultiplier tube, at the expense of the response time at high load values. The system allowed visualization of individual pulses as well as summation of pulses given in sequence. Using a virtual oscilloscope, we were able to record and sum the integral over as many successive shots as desired.

CONCLUSIONS

Time-resolved dosimetry of radiation requires fast responding electronics and sensors, pulse-to-pulse recording and real-time visualization of transients signals using ≥ 200 MHz bandwith digital oscilloscopes. Once this prerequisite is fulfilled, the machine operator committed to an experiment involving FLASH irradiation has to anticipate which dose the LINAC is going to deliver to the target surface. This is done by programming the beam parameters including the energy, peak current, pulse width, repeat frequency and pulse number which altogether determine the beam-on time and the integral dose at a given distance from the source. In our experience a repeatibility of ± 1.2% can be reached in routine with a LINAC in good condition and carefully maintained at thermal equilibrium.

A control of the incident dose should also be made a posteriori. It requires techniques that, with the exception of real-time integration of the electron charge collected at the primary graphite tore, are based on optical detection. We explored the linearity of these techniques in a range of peak dose-rates from 2x10 2 to 4x10 7 Gy.s -1 . These techniques include dosimetry in water based on formation of the methyl viologen radical cation, and Cerenkov light emission. The methyl viologen radical allowed us to analyze the response of Gafchromic™ EBT3 films in terms of dose-rate. These films are widely used in radiotherapy for control of the dose and field shape. We provide here direct experimental evidence that EBT3 films give the same response over a wide range of dose-rate, from 1 mGy to over 30 Gy per microsecond pulse with ±2% precision. EBT3 films may therefore be used confidently for measurement of the entrance and exit dose applied to laboratory animals, or to determine the depth-dose profile in solid water stacks at dose-rates exceeding by several orders of magnitude the specifications of ionisation chambers. Solid scintillation was also explored with

the Lynx system at peak dose-rates of 2.4-3.5 10 6 Gy.s -1 to help in determination of the shape and homogeneity of large radiation fields.

Cerenkov light has been used by others authors as a variable-delay probe beam with picosecond time resolution [START_REF] Grigoryants | Pulse radiolysis system with picosecond time resolution referred to Cherenkov radiation[END_REF]. Here we show that the small size of the optical probe we made for time-resolved analysis of the Cerenkov light allows high spatial resolution and is useful for small-field dosimetry at both low and very high dose-rate. Our results show that this technique may substitute for ionisation chambers at low dose-per-pulse and to Faraday cups at high dose-per-pulse through a dynamic range ≥ 4 10 4 , namely, from 1 mGy up to 40

Gy in a single microsecond pulse. Therefore, this technique is appealing in the prospect of future machines for FLASH radiotherapy in the clinics. Interestingly, the Cerenkov light has recently been considered for the dosimetry of X-ray, electron and proton beams [START_REF] Teymurazian | Megavoltage X-ray imaging based on Cerenkov effect: a new application of optical fibres to radiation therapy[END_REF][START_REF] Glaser | Video-rate optical dosimetry and dynamic visualization of IMRT and VMAT treatment plans in water using Cherenkov radiation[END_REF][START_REF] Zlateva | Cherenkov emission dosimetry for electron beam radiotherapy: a Monte Carlo feasibility study of absolute dose prediction[END_REF].

Other systems were tested including MOSFETS and diamond dosimeters. MOSFETS failed in our application due to both a large drift with the temperature rise accompanying big pulses and a reduced lifetime beyond 100 Gy cumulative dose. Diamond dosimeters were subject to large dose-and dose-rate effects. None was able to provide a time-resolved profile of the electron pulses at high dose-rate.

Yet the techniques analyzed in this study warrant further developments to be accepted as standard dosimetric methods in clinical applications, our conclusion at the end of the day is that they give a proportional response with the dose over a wide dynamic range with an overall precision better than ± 2%. This notwithstanding, as shown here FLASH is subject to the same pitfalls as conventional dose-rate radiotherapy relating to small-field effects and 8.

  . Beam collimation was achieved by means of a 22-mm thick graphite tore (d = 1.76, Schunk Electrographite, Lauffen, Germany) bolted onto a high density polyethylene insulator. Various collimators with 10-50 mm circular bore drilled in the graphite were used. In experiments involving local irradiation of mice, a secondary graphite collimator with suitable aperture delimitating a precise radiation field, was set to protect animals from whole-body irradiation. At high beam current (≥ 10 mA, full emission lobe) time-resolved measurement of the beam current was carried out by grounding through a load resistor (R) the charges collected in the core of the primary graphite collimator. The signal was sent through a coaxial wire in parallel to a National Instruments PXI-1033 digitizer with PXI-5152 (300 MHz, 8-bits) and PXI-5105 (60 MHz, 12-bits) amplifiers and a Nicolet Instruments (Madison, Wisconsin) Pro 92 oscilloscope with 200 MHz, 8-bits and 20 MHz, 12-bits amplifiers. Both NI and Nicolet 8bits ADCs were set in a dual base-time, slave mode triggering the 12-bits amplifiers to receive and store signals from the LINAC as well as from the optical stage.

Figure 1 .

 1 Figure 1. General design and main components of the Kinetron linear electron accelerator. RF: radio-frequency. WG: wave guide. RFW: radio-frequency window. AFC: automatic frequency control. IPC: incident power coupler. RPC: reflected power coupler. W: electron exit window. TEG, triode electron gun. GC: graphite collimator. NI PXI-5105, NI PXI-5152, Nicolet: high-speed oscilloscopes. The dotted and continuous lines represent command and power lines, respectively.

Figure 2 . 1 :

 21 Figure 2. Setup for irradiation of samples. 1: Accelerator framework. 2: Electron beam output. 3: High density polyethylene insulator. 4: 22 mm-thick primary graphite collimator. 5: Secondary graphite collimator (in vivo studies). 6: Bench. 7: Ball-bearing sliding aluminum plate. 8: Coaxial wire to remote oscilloscopes. R: Leak resistor. The optical system for real-time detection of transient absorption in chemical dosimeters is shown in Figure 3.

Figure 3 . 9 :

 39 Figure 3. Optical system for time-resolved analysis of transmittance changes in liquid samples. 9: Pulsed 150-W Xenon or continuous 150-W Xenon-Mercury lamp. 10: 150-W Tungstenhalogen bulb. 11: Fused silica lenses. 12: Remotely controlled, total internal reflexion prism. 13: Solenoid-activated light shutter and long-pass edge filters. 14: slit, cuvette holder and fused silica, gas-tight cuvette (10-or 20-mm path length). 15: UV-enhanced plane mirrors. 16: Holographic grating monochromator. 17: Fast response photomultiplier. 18: High-voltage supply. 19: Load resistors. 20: Back-off amplifier. 21: Digital oscilloscope. The figure is not to scale; the total light path is 5.9 m.

Figure 4 .

 4 Figure 4. In-depth calibration of the Markus ionization chamber and light absorption spectrum of Gafchromic™ EBT3 films. A. Variation of the charges collected at the Advanced Markus ionisation chamber irradiated in the CONV mode as a function of the depth in a PMMA stack. B. Light absorption spectrum of irradiated Gafchromic™ EBT3 films at 24 h post-irradiation.

Figure 5 . 1 :

 51 Figure 5. Exploded view of the probe used for collection and dosimetry of Cerenkov light. 1: duraluminum tube. 2: adaptor lens and washers. 3: fused silica cylinder. 4: aluminum foil. 5: hollow cap. 6: optical fiber bundle.

Figure 6 .Figure 7 .

 67 Figure 6. Screenshots and computer analysis of the electron current collected at the primary graphite collimator.We used a graphite collimator with a 40 mm wide bore and a load resistor of 50 Ω. Except in D, the peak current (full emission beam) was set at 103 mA, 30% of which was collected by the primary collimator. A: Oscilloscope time profiles of single, isolated pulses of width comprised between 0.10 and 2.20 µs. The area under the peak (in units of Vxs) for any given pulse allows determination of the electron fluence. B: Linearity of the fluence with the pulse width (single pulses). Spearman's test: y = ax (a = 1.542 ± 0.002 V), R = 0.999, p < 10 -5 . C: Linearity of the fluence with the pulse number in an experiment in which pulses (1.00 µs pulse width) were fired at 100 Hz repeat frequency. Spearman's test: y = ax (a = 1.583 ± 0.001 µV.s), R = 0.999, p < 10 -5 . The signal integrals (B, C) correspond to the summation

Figure 8 -

 8 Figure 8 -Small-field effects behind a secondary collimator. Left, variation of the dose with the distance between the target and the rear face of the secondary collimator. Right, dosimetric profile of the irradiation field along the horizontal (H) and vertical axis (V) as a function of the distance (mm, color code at upper left) to the collimator. Details are given under § 3.1.3.

Figure 9 .

 9 Figure 9. Use of the MV •+ radical for dosimetry at FLASH dose-rate. A. Formation of the MV •+ radical in N2O-saturated formate buffer following a 2 µs pulse of 15.5 (bottom) or 30.0 Gy (top). The optical signal was received at the monochromatorphotomultiplier stage during a pulse of the Xenon lamp synchronized with the electron pulse. Each experiment was repeated once to demonstrate the reproducibility of the process. B. MV •+ radical-based dosimetry of pulses in sequence. A 3 mM solution of MV 2+ in N2Osaturated formate buffer was exposed to seven 0.25 µs pulses of electrons at a repeat frequency of 100 Hz with online optical detection at 603 nm. The integral dose delivered to the irradiated sample was 26.1 ± 0.1 Gy. C. MV •+ radical-based demonstration of the inversesquare law fitting the intensity of the radiation as a function of the distance (r, mm) to the electron output window. The accelerator was set at 150 mA, 1.0 µs pulse width. Spearman's test: y = ax (a = 0.692 ± 0.012 Gy.m 2 ), R = 1.000, p < 10 -7 .

Figure 10 .

 10 Figure 10. Time-course of the decomposition of the MV •+ radical. A. The MV •+ radical formed during exposure to a 2.0 µs pulse of 31.2 Gy decayed slowly in the course of a biphasic process. The data points (two separate experiments) were fitted to the biexponential time-dependent equation (1) (see § 3.2). B. Effect of the instability of the methyl viologen radical under protracted irradiation at a low dose per pulse. We used equation (2) to calculate the evolution of the optical density of a methyl viologen dosimeter exposed to N pulses, 1.5 mGy each at a repeat frequency of 200 or 30 Hz. The black line represents the theoretical absorbance A0 = G * D * ε that would be obtained if the MV •+ radical

dose distribution in the depth of irradiated tissues. These heterogeneities can be modeled on a case-by-case basis with the aid of Monte-Carlo simulation [5].

ACKNOWLEDGMENTS

This work is dedicated to the memory of Alexis Sacharidis who designed hardware integration of dosimetric devices with NI LabView™.

FUNDING

This work was supported by the Agence Nationale de la Recherche (ANR-14-CE36-0008 program), the Comprehensive Cancer Center "SIRIC" program (INCa 2011-189) of Institut Curie and the Nanotherad IDEX (Paris-Saclay University). Generous financial aid from Mrs.

Marie Mahler (Gisikon, Switzerland), is also acknowledged.

was stable. The drop in aborbance relative to A0 after 15,000 pulses delivered at 200 Hz (T = 75 s) or 30 Hz (T = 500 s) amounts to 8.2% and 21.3%, respectively. Gy.s -1 ) and CONV modes (≤ 0.01 Gy.s -1 ). Dosimetry at both dose-rates was performed as described in text ( § 3.3). Data were blinded to mask any effect from the tester's expectation.

The optical densities measured at 634 and 583 nm were drawn for best fit to a biexponential equation, Table 1 Main characteristics of the Kinetron.

Output energy (MeV) 1 Pulse repetition frequency (Hz)

Pulse width (µs) Peak beam current (mA) 2 Maximum dose in a 1.0 µs pulse at 1 m (Gy) 3 Minimum dose per pulse at 1 m (Gy) 3 Mean linear energy transfer (LET, keV/µm) Relative biological effectiveness (RBE) 3.7 -5.1 1 -250 0.05 -2.20 0.04 -250 1.7 0.001 0.19 1.00 1 According to the Joule's law, the energy of output electrons decreases as the peak current during the pulse increases. For animal studies the kilovoltage applied to the magnetron (DQing) was adjusted to deliver a constant output energy (4.5 ± 0.1 MeV) in order to provide the same in-depth dose profile whatever the peak current used.


S Plaszczynski

J Peloton

C Arnault

J E Campagne

Analysing billion-objects catalogue interactively: Apache Spark for physicists

Keywords: Large-scale structure of Universe, galaxies:statistics, catalogues, distributed programming languages

Apache Spark is a Big Data framework for working on large distributed datasets. Although widely used in the industry, it remains rather limited in the academic community or often restricted to software engineers. The goal of this paper is to show with practical uses-cases that the technology is mature enough to be used without excessive programming skills by astronomers or cosmologists in order to perform standard analyses over large datasets, as those originating from future galaxy surveys. To demonstrate it, we start from a realistic simulation corresponding to 10 years of LSST data taking (6 billions of galaxies). Then, we design, optimize and benchmark a set of Spark python algorithms in order to perform standard operations as adding photometric redshift errors, measuring the selection function or computing power spectra over tomographic bins. Most of the commands execute on the full 110 GB dataset within tens of seconds and can therefore be performed interactively in order to design full-scale cosmological analyses.

Introduction

In 2002 Google released the mapReduce framework (see e.g, [START_REF] Dean | MapReduce: simplified data processing on large clusters[END_REF] a new parallel programming model exploiting efficiently many data centre hardware. In 2006 its open source implementation and many other tools related to the growing field of "Big data" emerged within the Hadoop ecosystem 1 .

In 2009 a research project started at UC. Berkeley [START_REF] Zaharia | Resilient distributed datasets: A fault-tolerant abstraction for in-memory cluster computing[END_REF][START_REF] Zaharia | Spark: Cluster computing with working sets[END_REF] to overcome some limitations met by Hadoop. In the next years, it captured many companies attention due to order of magnitude better performances on huge distributed data sets. It is today a very famous and active open source project named Spark owned by the Apache foundation2 and used by more than 1000 companies worldwide.

Spark was developed essentially in scala, a multi-paradigm language that appeared in 2004 within the java ecosystem, that revived functional programming (FP) 3 . FP deals with some old concepts related to λcalculus where functions become the central parts of programming. This approach allows the lazy evaluation mechanism which will be discussed in Sect.2. No knowledge of FP or scala is necessary to use Spark. Its API exposes efficient bindings to the python or R languages.

While under the influence of business and web companies, Big Data technologies were emerging, the data trend in astronomy was also rising rapidly. Typical surveys today reach Ter-abytes of raw data and Petabytes in a very near future (e.g., [START_REF] Zhang | Astronomy in the big data era[END_REF]. On the cosmological side, even after the reduction of images to user catalogues the BOSS experiment4 already imaged millions of galaxies. While the next generation of spectroscopic survey DESI5 is planned to image tens of millions of galaxies, the DES photometric survey6 already provides such a data volume. In the next decade LSST7 should reveal the properties of billions of galaxies that will be collected over 10 years. Moreover the study of the galaxy distribution properties requires tens to hundred times more simulated samples (mock catalogs).

We have also entered the Big Data area in astronomy not only for image treatment but also for analysing the output catalogues, which will be more and more limited by I/O throughput. Fortunately, such amounts of data are not impressive for today Spark industry standards.

Some exploratory work on using Spark in astronomy has already been shown by experts who ingested some large astronomical catalogues and organized them in a way where very optimized queries (as catalogues cross-match) can be performed efficiently [START_REF] Zečević | Introducing AXS: A framework for large-scale analysis of astronomical data[END_REF]Brahem et al., 2018a;[START_REF] Carretero | CosmoHub and SciPIC: Massive cosmological data analysis, distribution and generation using a Big Data platform[END_REF][START_REF] Riccio | C 3 , A Command-line Catalog Cross-match Tool for Large Astrophysical Catalogs[END_REF]. Although very promising for the future, this database-driven approach shows only one aspect of Spark possibilities.

Indeed Spark is a large framework and Fig. 1 shows its core modules, which can be complemented by numerous li-Apache Spark SQL Streaming GRAPH X

Machine Learning

Figure 1: Apache Spark provides four different modules. The SQL module focuses on the manipulating structured data such as tables. The Machine Learning module provides tools to perform distributed machine learning. The Streaming module allows to write streaming application. The GraphX module is the part of the API for graphs and graph-parallel computation. In this paper, we focus only on the Apache Spark SQL module.

braries. For this paper we will focus essentially on the SQL part, allowing users to perform efficient queries and analyses on catalogs. But Spark also offers machine-learning possibilities (MLLIB module [START_REF] Meng | MLlib: Machine Learning in Apache Spark[END_REF]) which can be complemented by external libraries 8 . A first application of the graph analysis module (Graph-X) to investigate the topology of cosmological simulations was recently presented in [START_REF] Hong | Constraining Cosmology with Big Data Statistics of Cosmological Graphs[END_REF], while Streaming possibilities open exciting perspectives for the real-time treatment of alerts by brokers in astronomy 9 . An introduction on using these different modules can be found in [START_REF] Zečević | Spark in Action[END_REF].

Apache Spark was also used to perform distributed processing of astronomical images (Zhang et al., 2015;[START_REF] Wiley | Astronomical Image Processing with Hadoop[END_REF] or to produce mock catalogues from huge N-body simulations [START_REF] Carretero | CosmoHub and SciPIC: Massive cosmological data analysis, distribution and generation using a Big Data platform[END_REF]. This illustrates the potential of the framework when used for fundamental sciences.

Our goal here is not to report the impressive results obtained by Spark experts, but to address the needs of end-users (physicists) who have little or no knowledge about Spark. We will not elaborate on the theory that can be found in many articles or textbooks (a classical introduction being [START_REF] Karau | Learning Spark. O'Reilly[END_REF] but try to show with concrete examples how this framework offers a simple way to achieve excellent results for analysing today's and upcoming catalogs in astronomy on computing clusters. Today's main language for astronomical data analysis being python, we will focus on that language. We will use some familiar concepts as dataframes and emphasise how external python codes can be reused. The use-case we develop is the following: one obtains a (large) set of astronomical files, not 8 see for instance https://maxpumperla.com/elephas. 9 as in https://fink-broker.readthedocs.io/en/latest/ necessarily produced by the user, and would like to have a look at some of the variables and perform some simple analysis over the entire dataset. As an example, we choose to use a catalogue of 6 billion of galaxies generated with a fast simulation and show how to perform easily an interactive tomographic analysis in a very reasonable time. Since a very common format in astronomy is FITS, we first developed an efficient connector named spark-fits10 described in [START_REF] Peloton | FITS Data Source for Apache Spark[END_REF], hereafter called SparkFITS18.

After presenting in a physicist's language what Spark is and what its advantages are for data analysis in Sect. 2, we describe how to start using Spark in python in Sect. 3. Then we design and optimize several Spark commands in Sect. 4 with increasing complexity leading to a full tomographic analysis over a simulated set of 10-years of LSST-like galactic data (Sect. 4.3.7). Several options are explored and performances discussed. AppendixA presents a more realistic treatment of photometric uncertainties.

What is Spark and why to use it

The main trend for analysing large astronomical datasets in the last two decades was to use High Performance Computing (HPC) approaches on super-computers. The challenge in using such a framework consists in optimising the arithmetic efficiency (number of data moves over number of operations) on more and more complex architectures and requires a high degree of skills. This can be performed with the Message-Passing-Interface (MPI, [START_REF] Gropp | A high-performance, portable implementation of the mpi message passing interface standard[END_REF]) but requires a careful tuning of the code for data exchanges between the nodes.

Starting with mapReduce, a new approach to address these problems has focused on High Through-put Computing (HTC) in data centres. Accent is put interacting with nodes that are "close" to the data and minimizing the traffic between them.

A key requirement for this is robustness and reliability, which is now achieved by performing computation in a fault-tolerant way on hardware which may not be fault-tolerant. A more in-depth overview can be found in [START_REF] Stickley | NebulOS: A Big Data Framework for Astrophysics[END_REF].

Spark is now for many years the leading solution in this domain. Hosted by the Apache foundation, the code is opensource and its community very active. In very general terms, "Apache Spark is a cluster computing platform designed to be fast and general" [START_REF] Karau | Learning Spark. O'Reilly[END_REF]. Although available for personal usage on laptops, it will reveal its full potential when running on a cluster on some large amount of data11 within a distributed file system.

To introduce the necessary vocabulary, we outline some of the features of the Spark solution and refer the reader to [START_REF] Karau | Learning Spark. O'Reilly[END_REF]; [START_REF] Tang | A Survey on Spark Ecosystem for Big Data Processing[END_REF]; [START_REF] Peloton | FITS Data Source for Apache Spark[END_REF] for a more thorough technical introduction.

In Spark, the user program, called the driver, launches multiple workers. Each worker can hold several executors (processes launched for an application) that execute a number of tasks.12 Workers read data by blocks (possibly from a distributed file system) and can persist computed objects in memory. In a typical work flow the Spark scheduler sends a set of tasks to the workers in order to compute locally over the nodes holding the data or in memory if they have been cached. The reduced outputs are then returned to the driver that can perform some final combination.

Although quite an evolved process, on a practical level the user only needs to interact with RDDs (Resilient Distributed Datasets, [START_REF] Zaharia | Resilient distributed datasets: A fault-tolerant abstraction for in-memory cluster computing[END_REF]) which are abstractions representing a distributed collection of objects. They provide a low level access to Spark features. Scientists will typically prefer to work with DataFrames (from the SQL module), which are built on top of RDD but with data organized into named columns and which also include the benefits of the Spark SQL's execution engine (Sect. 4.3.1).

Using Spark presents a number of advantages that are particularly well suited to data analysis tasks in the physical sciences, as we discuss below.

Simple parallelisation. The use of the mapReduce paradigm allows a coarse level of parallelisation without ever writing complicated code or directives. This is best illustrated by an example. Suppose you are looking for the maximum value of a huge set of numbers. The driver requires workers to compute their maximum value. The results from this map step are then returned and combined in the reduce step, here by taking the maximum of their values. Such a behaviour could also be implemented with other mechanisms such as MPI, but would require writing some extra piece of code while, quite remarkably, it is builtin feature in Spark.

Optimization. In imperative languages like for example C, C++ or Fortran, each instruction generates some low level commands that are immediately executed. Spark uses the mechanism of lazy evaluation, meaning that most instructions (called transformations) do not immediately trigger execution, but instead updates the construction of a Direct Acyclic Graph or DAG which in physicists language represents a pipeline. It is often a common surprise to newcomers to see how fast some Spark commands seem to perform, but this is only because the pipeline has not yet been closed by specifying an action, so nothing was actually run (see an example in Sect. 4.3.2). Then, the python front-end to Spark can be seen as a way of declaratively specifying a calculation that is then passed to the Spark layer for evaluation. An interesting aspect of this is that knowledge of the full pipeline structure allows Spark to perform some precise optimizations based on a set of logical/physical rules before running it. As long as one makes use of Spark functions, the user is relieved of the burden of code optimisation such as accessing different parts of data in different places.

Cache and interactive analysis. Maybe the most interesting feature for the end-user analysis is the ability of putting the data in cache. While quite a complex operation behind the hood, on the user side this is as simple as calling the cache() function. Even when the cluster does not allow all the data to be held in memory, putting some part in cache and spilling the rest over disk can still be interesting (see Fig. 3 of SparkFITS18). A key point is that in general, you do not need to put all the variables in cache, only those that will be accessed several times. The memory requirements are then drastically reduced.

Once the relevant variables are cached, all computations can be carried out efficiently enough that an interactive analysis becomes generally possible. For demonstration purposes in the following sections, we have put 110 GB of data in memory on a small cluster, but one can achieve easily 1 TB on large centres like the National Energy Research Scientific Computing Center (NERSC) 13 .

Performances and scaling. Spark, as a largely used industrial standard, gives generally excellent performances for analysing large datasets. We do not claim that Spark will always give the ultimate performance as depending on the specific domain some implementations with MPI or using relational databases, carefully tuned by experts, may outperform it. But having at disposal such a simple tool to reach such high level of efficiency is a great opportunity for physicists facing Big Data challenges.

The attractiveness of Spark also resides in the fact that one can design some analysis of a partial set of data on a personal laptop and later port it to large data centres knowing that in most cases performance will scale.

Working with Spark

As we saw previously, Spark's basic abstraction is the RDD which represents a coarse-grained distributed collection of objects which allow in-memory fault-tolerant computations. Dataframes were later introduced in 2015 within the SQL module [START_REF] Armbrust | Spark sql: Relational data processing in spark[END_REF]. They add to the RDD the knowledge of its data structure which allows a higher level of optimization. Although one can retrieve an RDD from a dataframe, one should always try to use the latter in the analyses, especially for python applications. Indeed the RDD performances are much worse in python than in the native scala language. This is not any more the case when working with dataframes, for which similar performances are obtained in most cases as shown in Sect. 4.3.8.

Most of astronomers use today the python language for their analyses and we emphasize that pyspark can be run within ipython or a jupyter notebook. We will then essentially discuss the python interface although some comparison with scala performances will be later shown. In our opinion, the scala language remains an interesting option to be looked at and we highlight some of the reasons why it may be useful to astronomers below:

• there exist some Spark kernels in scala to work in jupyter notebooks14 ,

• the HEALPix package15 can be build in java (and therefore used in scala) although not including all functionalities,

• the JEP package16 allows to make calls to external python modules as numpy or matplotlib, and the library Py4J17 provides an efficient bridge between python and scala.

Concerning data formats, spark-fits allows the reading of FITS distributed binary tables and also in recent versions of images. Some connectors also exist to read the HDF5 format (see references in SparkFITS18).

A full interactive analysis

We present a suite of rather standard operations18 an enduser can perform using the output of a 10 years LSST simulation presented in Sect. 4.1. Commands were run interactively in the pyspark shell on a cluster described in Sect. 4.2. They are explicitly shown as boxed inline code with their output following right after. We benchmark different options in some cases and focus more on performances in Sect. 4.3.8.

Simulation

In order to work with physics-oriented data, we built a catalogue of galaxies using the CoLoRe fast simulation19 corresponding to 10 years of LSST data-taking. Point-like galaxies are generated in the z ∈ [0, 2.5] redshift range assuming a standard ΛCDM cosmology, with a selection function coming from the LSST/DESC 2pt validation working group20 and shown on Fig. 2.

The generated catalogue consists of 6 billions of galaxies with their types, RA/DEC positions, cosmological redshifts and the redshift-space distortion (RSD) displacements. The CoLoRe simulation was run at NERSC and produced 32 FITS files of 112 GB in total that were imported to our HDFS cluster [START_REF] Shvachko | The Hadoop distributed file system[END_REF]. By importing, we mean simply copied to the HDFS cluster without any post-processing. Then, spark-fits gives transparent access to all the Spark advantages.

This procedure is mainly pedagogical since it introduces latencies for transferring data to another cluster.

Whenever possible, one should instead attempt to produce and analyse the data on the same cluster, which would have been possible on the NERSC supercomputer since it offers a Spark API (more in SparkFITS18).

Infrastructure

The cluster we used for this work, located at Université Paris-Sud in France21 is rather modest in order to illustrate the fact that you do not need huge resources to achieve spectacular results. It consists of nine 36 GB machines, each with 18 cores, running over HDFS. Since the total memory fraction dedicated to the cache is set to 0.6, the usable memory amount ( 170 GB) is largely sufficient to hold our full dataset. Then, we used the following set-up

• 1 driver (4 GB RAM)

• 8 executors (ie. workers) each using 17 cores and 30 GB of RAM.

The python interactive shell is run as followed pyspark --driver-memory 4g \ --total-executor-cores 136 \ --executor-cores 17 \ --executor-memory 30g 4.3. Spark analysis 4.3.1. Using dataframes Dataframes in Spark come from the Spark.SQL module. Similar to the pandas one 22 , they can be viewed as a set of named columns over which you can perform operations, but in a distributed environment. Some native Spark functions act on them and should be used as much as possible since they have been very optimized. They are available through the following import command from pyspark . sql import functions as F Once substantial data reduction has been achieved, we can recover a standard pandas dataframe with the toPandas() method, which opens the door to further standard python analysis. It is worth mentioning that toPandas() returns to the driver memory, so one should keep in mind that the data reduction (reduce methods) must be performed before.

Reading the data

We begin by loading the (set of) FITS files using spark-fits. CoLoRe FITS format stores separately the cosmological ( Z_COSMO) and RSD (DZ_RSD) redshifts, but since we only want to work on their sum, we construct the z column on the fly: gal = spark . read . format ( " fits " ) \ . option ( " hdu " ,1 ) \ . load ( " hdfs : path / to / fits / directory " ) \ . select ( " RA " ," Dec " , \ ( F . col ( " Z_COSMO " ) + F . col ( " DZ_RSD " ) ) \ . alias ( " z " ) )

This represents a (Spark) dataframe object. We show two ways of accessing columns in the select function, either simply through their names (strings) or, when some operations are to be performed, through F.col(..) that returns Columns objects.

Other dataframe standard ways to access columns are through gal.RA or gal['RA']. Note the use of the "z" alias to rename the new column.

You can now print the dataframe schema: Following the principle of lazy evaluation, it is important to realize that, at this level, data is not (yet) physically loaded: only the FITS header is read and the DAG updated.

Adding photometric smearing

We would like to add now the effect of the photo-z (PZ) resolution of the instrument. For LSST this will be close to a Gaussian and the upper requirement for Large Scale Structure analysis is given in The LSST Dark Energy Science Collaboration et al. ( 2018) as σ z = 0.03(1 + z) . A lot of work based on template-fitting or Machine Learning techniques shows that the actual distribution is more complicated than a simple Gaussian. However to not distract the reader from the main goal of this paper, we have relegated to AppendixA a more realistic discussion of this effect and how it can be implemented within Spark.

For the Gaussian case, we simply add a column to the gal dataframe according to from pyspark . sql . functions import randn gal = gal . withColumn ( " zrec " , gal . z + 0 . 03 * ( 1 + gal . z ) * randn () )

. astype ( ' float ') )

Again, when this last command is executed the DAG is further filled but no data has been physically generated.

Let us now show 5 samples, which triggers an action:

gal . show ( 5) --------+---------+---------+--------- --------+---------+---------+---------+ only showing top 5 rows This happens within seconds (see Table 1). How is that possible? This is the very idea of lazy evaluation: if you only want to look at a few samples is it worth loading all (110 GB) of the data? Here Spark analyses the full pipeline, optimizes it and only physically read the very first block

+---------+---------+---------+---------+ | RA| Dec| z| zrec| +-

Caching

Now that we have defined which data we want to use in our analysis, we put them in cache. This is achieved with the cache() function. Some finer level of details can be obtained with the persist(level) function that allows to specify the storage level. cache() corresponds to persist(MEMORY_ONLY). You may use persist(MEMORY_AND_DISK) if your cluster does not have enough total memory. It was shown in SparkFITS18 that good performances can still be obtained in this case. Note that serializing the objects might also improve performances (e.g. level=MEMORY_ONLY_SER) but it was not observed in our case. In order to trigger caching, one must call an action as counting the total number of galaxies which requires access to the full data: print ( gal . cache . count () )

5926764680

In our case, putting all the data in cache by counting them takes about 90 s, with about 20 s coming from the PZ computation.

Getting some basic statistical information

Some basic statistical informations may be obtained on some (or all) variables with: +-------+-------------------+------------------ ------+-------------------+------------------+ If we only need some specific values (e.g minimum/maximum), it is more efficient to use the Spark functions: minmax = gal . select ( F . min ( " z " ) ,F . max ( " z " ) ) \

+ |summary| z| zrec| +-------+-------------------+------------------+ | count| 5926764680| 5926764680| | mean| 0.875229444425171|0.8752293689731887| | stddev|0.47360539092073933|0.4771461812884577| | min| -5.93947E-4| -0.12403674| | max| 2.4352543| 2.943411| +-
. first () zmin = minmax [ 0 ] zmax = minmax [ 1 ] Table 1 gives the measured user-time in each case: once the data are in cache, one observes that all those commands run in seconds.

Histograms

We now wish to go further and study the redshift distribution of galaxies. Although RDDs provides some command to build histograms, we will show that it is much more efficient to design a function using dataframes capabilities, as follows 1. add a new column to the dataframe containing the bin number, 2. group the data by this number, 3. count the number of values in each group, 4. sort the bin index by ascending order.

Adding the z bin number column (labelled "bin") is done most efficiently using standard column operations: Nbins = 100 dz = ( zmax -zmin ) / Nbins zbin = gal . select ( gal .z ,\ (( gal . z --dz / 2 ) / dz ) . astype ( ' int ') \ . alias ( ' bin ') )

Then, grouping by the bin column, counting its members and sorting in ascending order is performed by: h = zbin . groupBy ( " bin " ) \ . count () \ . orderBy ( F . asc ( " bin " ) )

Finally we may want to add the bin locations, drop the bin number and go back to the python world by recovering a standard pandas object: pd = h . select ( " bin " ,\ ( zmin + dz / 2 + h . bin * dz ) . alias ( ' zbin ') ,\ " count " ) \ . drop ( " bin " ) \ . toPandas ()

The histogram is obtained in about 10 s (Table 1) which is impressive for running on 6 10 9 data records.

We can now study for instance how the selection function varies with the PZ smearing which is shown on Fig. 3.

In order to prepare for the next part, let us see how to build the histogram by calling an external function. Operations may be applied onto dataframes with a User Defined Function (UDF) as: But we find that there are performance issues since execution time goes from the previous 10 s to 2 minutes.

To alleviate this issue, Spark introduced recently (v2.3.0) pandas udf's. Previously, standard UDF's were processing one row-at-a-time, which turns out to be slow for large numbers of rows. Newly introduced pandas udf's are built on top of Apache Arrow23 which, among several other features, bring vectorized optimization. Scalar pandas udf's24 are used for vectorizing scalar operations using pandas.Series (groups of rows). The UDF is re-written as: import pandas as pd from pyspark . sql . functions \ import pandas_udf , PandasUDFType @pandas_udf ( " float " , PandasUDFType . SCALAR ) def binNumber ( z ) : return pd . Series (( z -zmin ) / dz ) zbin = gal . select ( gal .z ,\ binNumber ( " z " ) . astype ( ' int ') \ . alias ( ' bin ') )

The user-time becomes 40 s which is better although not optimal and will be discussed more in Sect. 4.3.8. The main lesson from this part for python users is to always work with dataframes and whenever possible with the native Spark.SQL functions.

Tomography

Measuring galactic over-density power-spectra over some redshift bins (called tomographic "shells") is a nearly optimal method in cosmology to study galaxy clustering, especially for photometric surveys [START_REF] Crocce | Modelling the angular correlation function and its full covariance in photometric galaxy surveys[END_REF][START_REF] Asorey | Recovering 3D clustering information with angular correlations[END_REF]. Measuring the cross-correlation between nearby shells gives also access to Redshift-Space-Distorsions even in photometric surveys where the radial information is strongly suppressed [START_REF] Ross | Measuring redshift-space distortions using photometric surveys[END_REF]. Cross-correlation between far-away bins is also of interest: since, neglecting magnification lensing, no cosmological signal is expected there, any observed correlation singles-out some remaining systematics (as PZ distribution tails). Such studies can be efficiently performed with Spark. We have chosen 10 redshit bins marked out as vertical lines on Fig. 4. Since the selection acts on the observed space, each bin receives contribution from the true redshifts according to the coloured distributions.

For each tomographic shell, we use the well-known HEALPix sphere tessellation 25 , which partitions hierarchically the sphere (using the nside parameter) with equal-area and iso-latitude pixels. This latter property allows an efficient computation of power-spectra over the sphere. Using this scheme, one can map each (RA/DEC) position on the sky to a single pixel number. So, the construction of the projected galaxy number onto a HEALPix map (nside = 512) is as easy as building an histogram but using the external ang2pix function to determine the pixel number. For performances we use pandas udf to call this function.

For a z ∈ [z 1 , z 2 ] shell import pandas as pd import numpy as np import healpy as hp nside = 512

# define the UDF @pandas_udf ( ' int ' , PandasUDFType . SCALAR ) def Ang2Pix ( ra , dec ) : theta = np . radians ( 90 -dec ) phi = np . radians ( ra ) return pd . Series ( hp . ang2pix ( nside , theta , phi ) )

# build the shell shell = gal . filter ( gal [ ' zrec '] . between ( z1 , z2 ) )

# build the HEALPix map map = shell . select ( Ang2Pix ( " RA " ," Dec " ) \ We end up with a standard HEALPix map on which one can perform further analysis. An example is shown on Fig. 5. We note that the python packages must be available on each executor.

z [0.61, 0.82] 250 350 Concerning performances (Table 1), each shell projection is obtained in about 30 s, quite independently of the galaxy population. All the 10 shells are obtained in about 5 mins.

From the maps, one can then compute auto and cross spectra using standard healpy functions. We illustrate some results we obtain for a few of them on Fig 6.

Performances

We already discussed the user-time measured for each step. They are summarized in Table 1. But are we far from the best possible ones in Spark? The native Spark language is scala which generally leads to the best performances. So, we have recoded and run all the previous commands in scala and compare performances in Table 1.

The initialization phase (load(HDU)) is slightly longer in scala than in python. Caching the 6 10 9 data records which is the most demanding part but only needs to be performed once, is obtained in our case in about 1.5 mins in each case. Then, the statistics part is slightly more efficient in python. As long as native dataframe operations are used both implementations yield a similar 10 s value for histograms. The use of an external (UDF) function is severely penalized in python. Using pandas udf's the user-time is reduced to 40 s which becomes reasonable although a factor of 3 higher than the corresponding scala implementation. The same kind of factor is observed

C i × j 1e 5 0 × 9 1 × 8 2 × 7 3 × 6 4 × 5 4 × 4
Figure 6: Tomographic power spectra reconstructed by cross-correlating the over-density maps with the bin numbering shown in the legend (see also Fig. 4). As discussed in Sect. 4.3.7 we see no power for well separated bins, some small contribution from the 4 × 5 adjacent one and strong power for the 4 × 4 auto-spectrum. when building the shells. Still, reconstructing a tomographic map in about 30 s remains very satisfactory for interactive work.

Let us compare these performances to the ones we would have obtained in a more standard imperative way. For instance, let us consider how we would have computed the min/max values of the PZ redshifts:

1. each FITS files is opened and read, 2. redshift uncertainties are added shooting random numbers, 3. the min and max values are computed and stored, 4. all min/max values are combined to obtain the lowest (highest) min (max).

We implemented that process in python: it takes about 45 mins to process the 32 files. Which is to be compared to the 2.8 + 12.4 + 97.7 + 1.8 2 mins we obtained with Spark. This sequential comparison is unfair since one could also implement some MPI logic to read the files in parallel on each node, or even post-process the files to distribute them on each core. The point is not getting the very best performance but already an excellent one with very little effort. Furthermore, while in the imperative approach all the data in memory is lost when the program ends, using the cache in Spark offers the possibility to investigate several aspects of the data in a single session without ever reloading them.

Conclusion

In the current days of ever-growing astronomical data, it is worth investing in a technology which gives physicists access to the interactive analysis of billions of objects, boosting timeto-physics considerably.

We have illustrated in this pedagogical introduction to users, how a 6 billion-objects realistic galactic catalogue could be analysed interactively producing histograms in about 10 s and tomographic bins in 30 s. Our main point is about simplicity. The access to a large volume of data in a distributed environment requires today quite involved technical skills that are relieved by simply using Spark (python) functions.

We emphasize that Spark is a Big Data technology, meaning it makes no sense using it for data volumes below a few tens of Gigabytes as one would lose its main computational advantages by paying its overhead costs. Although it is convenient to develop the software on a personal computer do not expect interesting performances there. The Spark technology is intended to be run on data-centres where the storage is a crucial part. However one does not necessarily need a very large cluster or a super-computer for that; we achieved very satistactory results on 8 machines. We focused here on rather simple operations without ever worrying where data reside on the cluster nor how they are organized, a process known as data partitioning. For performance, one always needs to minimize network connections among the workers and this may need to be investigated in some cases, as for instance when cross-matching catalogs [START_REF] Zečević | Introducing AXS: A framework for large-scale analysis of astronomical data[END_REF]Brahem et al., 2018b). It can be however a heavy process since data needs to be shuffled among the different executors, and possibly re-written to the disk. It further depends on the type of query that is to be optimized. Spatial locality may be achieved straightforwardly in 2D with a Healpix indexing and in 3D with an octree26 . But some more evolved indexing may be required for instance when also considering the color or shape of galaxies. This is an important topic, but one that goes beyond the scope of this paper and will be discussed in a forthcoming one.

We have addressed the question of interfacing pyspark to some external python code, which can be done reasonably well using (pandas) user-defined-functions. This however requires I/O's to be quite simple. A difficulty we encountered was that building some application always requires access to some external python libraries that must be available on each executor. While Spark offers some integrated mechanism to transmit binary modules to all workers, it is not a very efficient solution. The administrator may install some set of common scientific libraries on each node, but it is not a very flexible solution. An elegant solution we found for NERSC was to build a dedicated shifter image (a docker variant) including Spark and the necessary libraries. This is a simple solution that can be put in place by each user, but running a docker-like image over each node is not always supported on every data centre.

One can also interface some C, C++ or Fortran codes with scala and the JNA (Java Native Access) library. Some examples were reported in Zhang et al. (2015); [START_REF] Wiley | Astronomical Image Processing with Hadoop[END_REF] and the interested reader may start with https://github. com/astrolabsoftware/Interfaces for a practical usage. On super-computers one may take benefit from combining the HTC and HPC approaches, a road that has still to be explored. This is an exciting perspective, that would allow to add science to a technology born once in a public lab. It is the goal of the Astrolab27 organization to have everyone interested in such a project join it. confusion may arise when matching the catalogue which results in tails in the inferred redshift values. Knowing precisely the photometric redshift distribution is an important ingredient to several key cosmological analyses as those related to weaklensing [START_REF] Mandelbaum | Weak Lensing for Precision Cosmology[END_REF].

In order to assess a more realistic LSST use-case and show how to implement it within Spark, we use some recent results from [START_REF] Choyer | Impact of photometric redshifts on the BAO scale determination in the LSST survey[END_REF]. The authors have reconstructed the distributions of the observed photometric redshifts given the true ones P(z rec |z) applying a full template-fitting procedure. Fig . A.7 shows some of these distributions which reveals that their main part is more packed than the Gaussian model but exhibits complex tails.

Our aim is, for our simulation, to shoot some random number according to these distributions. We use for that a standard Monte Carlo technique: to generate samples according to some target distribution P, one generates some uniformly distributed samples and transform them according to the inverse of the P cumulative distribution, or, in statistical language u ∼ U(0, 1), F -1 (u) ∼ p when F(z) = z -∞ P(x)dx. The cumulative function can be computed and inverted numerically.

In our case, the authors of [START_REF] Choyer | Impact of photometric redshifts on the BAO scale determination in the LSST survey[END_REF] have provided a table of 300 photometric distributions according to their true redshifts in the z ∈ [0, 3] range. For each true value z i we computed the cumulative sum of the corresponding photometric distribution P(z rec |z i ) and numerically invert it. Since in the inversion one looses the abscissa regularity, we re-sampled the values onto a linear grid of 1000 points over the [0,1] interval. We then stored the resulting values in a table 300 × 1000 as a plain text file. Then, for each true redshift value z i and an associated uniformly generated number u i , one just need to read the proper value from this table. We have checked that the full procedure applied to the gaussian case reproduces the results shown in Sect 4.3.6.

The implementation of the method within Spark is now rather simple:

1. read the inverse-cumulative table 2. add to the dataframe a column of uniformly generated numbers ("u"). 3. for each true redshift ("z") and random sample ("u"), locate the index within the table and get the value from a user-defined-function (pandas udf).

The code reads as follows # add column of uniform random numbers gal = gal . withColumn ( " u " ,F . rand () )

# t r a n s f o r m with the inverse -c u m u l a t i v e table gal = gal . withColumn ( " zrec " , z_phot ( " z " ," u " ) + dz / 2 ) \ . drop ( " u " ) # do not need u anymore

Using the method discussed in Sect. 4.3.6, we histogram the "zrec" column which gives the selection function shown on Fig.

A.8. The non-Gaussian behaviour of the photometric distribution dramatically affects the selection function in a non-trivial way, which to our knowledge is shown for the first time.

Figure 2 :

 2 Figure 2: Density number of galaxies (black, left axis) and their cumulative sum (dashed-blue, rights axis) generated by CoLoRe corresponding to 10 years of LSST data-taking.

  RA: float (nullable = true) |--Dec: float (nullable = true) |--z: float (nullable = true)

  gal . describe ( [ 'z ' , ' zrec '] ) . show ()

Figure 3 :

 3 Figure 3: Difference of galactic density when applying a Gaussian photometric smearing.

Figure 4 :

 4 Figure 4: Position of the 10 tomographic bins in the observed space (vertical black lines with the numbering convention up) and contributions from the true redshifts to each bin assuming a Gaussian PZ smearing with σ z = 0.03(1 + z) (filled curves).

  25 https://healpix.jpl.nasa.gov/ . alias ( " ipix " ) ) \ . groupBy ( " ipix " ) \ . count () \ . toPandas ()#back to python world myMap = np . zeros ( 12 * nside ** 2 ) myMap [ map [ ' ipix '] . values ] = map [ ' count '] . values

Figure 5 :

 5 Figure 5: Mollweide projection of the contrast density map N- N N obtained from our LSST simulation on bin 4 (see Fig. 4).

Figure

  Figure A.7: Distributions of some photometric redshift distributions given the true values shown in the box, as obtained in Choyer et al. (2019) using a full template-fitting method. Note that we do not use a quality cut. The dashed lines corresponds to the simplified Gaussian model with σ z = 0.03(1 + z).

Figure A. 8 :

 8 Figure A.8: Histogram of the reconstructed redshift of galaxies in our 10-years LSST sample (full line) compared to the underlying real one (dashed), using the realistic photometric distributions from Choyer et al. (2019).

  import numpy as np import pandas as pd # read the inverse -c u m u l a t i v e file cuminv = np . loadtxt ( ' cum_inv . txt ') # we know the binning that were used dz = 0 . 01 du = 1 / 1000 . # find indices and return the table values @pandas_udf ( ' float ' , PandasUDFType . SCALAR ) def z_phot ( zr , u ) : iz = np . array ( zr / dz , dtype = ' int ') iu = np . array ( u / du , dtype = ' int ') return pd . Series ( cuminv [ iz , iu ] )

Table 1 :

 1 User-time (in seconds) for the various analysis steps described in the text using the python commands (first column) and the scala ones (second column). Results were obtained by running sequentially the commands from the top to the bottom of the table, 10 times and were then averaged.

	Section analysis	python	scala
	4.3.2	load(HDU) PZ + show(5)	2.8 ± 0.1 12.4 ± 0.6 13.7 ± 1.2 8.8 ± 0.2
	4.3.4	cache (count)	97.7 ± 4.0 95.4 ± 5.0
		stat(z)	3.9 ± 1.5	4.9 ± 2.5
	4.3.5	stat(all)	9.8 ± 1.0	11.0 ± 0.9
		minmax(z)	1.8 ± 0.3	3.2 ± 0.7
		histo (dataframe)	11.5 ± 1.5 13.0 ± 0.8
	4.3.6	histo (UDF)	114.9 ± 5.6 13.9 ± 1.2
		histo (pandas UDF) 43.3 ± 4.5	-
	4.3.7	1 shell	30 ± 3	13 ± 2
		all shells (10)	307 ± 34	130 ± 18

https://hadoop.apache.org/

http://spark.apache.org/

FP should not be compared to procedural or object-oriented programming which are both imperative languages.

http://www.sdss3.org/index.php

http://desi.lbl.gov

https://www.darkenergysurvey.org

https://www.lsst.org

https://astrolabsoftware.github.io/spark-fits

by large we mean beyond several tens of Gigabytes.

In the examples below, we restrict the case to one executor per worker, with several tasks.

https://www.nersc.gov

e.g. https://toree.incubator.apache.org

http://healpix.sourceforge.io

https://pypi.org/project/jep

https://www.py4j.org/

A jupyter notebook is available at https://github.com/ astrolabsoftware/1807.03078

https://github.com/damonge/CoLoRe

https://github.com/LSSTDESC/2pt_validation

https://www.informatique-scientifique.u-psud.fr/ services/spark.html

https://pandas.pydata.org/

https://arrow.apache.org/

There are two types of pandas udf's: scalar and grouped map. In this paper we focus only on the former.

https://astrolabsoftware.github.io/spark3D offers a simple dataframe interface for doing it

https://astrolabsoftware.github.io

Acknowledgements

We acknowledge the use of the HEALPix package (Górski et al., 2005) and the CoLoRe fast simulation with support from David Alonso, Anze Slosar and Javier Sanchez that we kindly thank. The Spark work was performed at the VirtualData center at Université Paris Sud and we thank Adrien Ramparison for the upgrades and maintenance of the cluster. The CoLoRe simulation was run at the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

//github.com/

AppendixA. Handling realistic photo-z distribution

We have seen in Sect 4.3.3 how to implement some very simple Gaussian smearing of the redshift variable. Beyond this simple approximation, the more complex situation is the following. In brief, a photometric instrument like LSST measures integrated fluxes in a few frequency bands (6 for LSST) and from these few values one infers the redshift. This is performed with the help of some external catalogues where the spectral energy distributions (SED) of a sample of galaxies has been precisely measured, for instance with a spectrometer. Since the effect of redshift is to shift the SED, the idea is to determine the redshift (noted z rec ) by measuring the displacement of the observations to match the catalogue SED. This can be performed with template-fitting or Machine Learning algorithms (see for instance [START_REF] Salvato | The many flavours of photometric redshifts[END_REF]; [START_REF] Pasquet | Photometric redshifts from SDSS images using a convolutional neural network[END_REF]). But, some