3D diffractive imaging of nanoparticle ensembles using an X-ray laser

We report the 3D structure determination of gold nanoparticles (AuNPs) by X-ray single particle imaging (SPI). Around 10 million diffraction patterns from gold nanoparticles were measured in less than 100 hours of beam time, more than 100 times the amount of data in any single prior SPI experiment, using the new capabilities of the European X-ray free electron laser which allow measurements of 1500 frames per second. A classification and structural sorting method was developed to disentangle the heterogeneity of the particles and to obtain a resolution of better than 3 nm. With these new experimental and analytical developments, we have entered a new era for the SPI method and the path towards close-to-atomic resolution imaging of biomolecules is apparent.

The determination of the structures of biomolecules at atomic resolution requires bright sources of radiation, which are unfortunately also energetic enough to degrade the object under observation [1]. All approaches to structure determination are primarily dedicated to overcoming, or working around, the effects of this radiation damage. In X-ray crystallography, large numbers of aligned molecules amplify the diffraction signal that can be obtained within the exposure that the sample can tolerate. The tolerable dose can be increased somewhat by cooling the crystals to cryogenic temperatures. Such cooling also allows electron microscopy-where the ratio of the image-forming to damage-causing radiation is more favourable-to record faint and noisy images of many uncrystallised molecules, which can then be used to build up a three-dimensional image. The extreme intensity and ultrashort pulses of X-ray free electron lasers (XFELs) potentially offer another way to obtain structural information from single macromolecules, but without the need for cooling [2]. Pulses of femtosecond duration can outrun radiation damage and essentially freeze the molecule in time [3,4].
Single particle imaging at XFELs consists of collecting coherent diffraction patterns from individual particles intersecting bright XFEL pulses. Theoretical work predicts that currently available XFEL sources generate enough scattered photons from single macromolecules to solve for their unknown orientations and reconstruct 3D structures of large reproducible biomolecules [5][6][7]. Proof-of-principle SPI experiments on biological particles [8][9][10][11][12][13][14] have highlighted the challenges of the approach i.e. the recording of a large number of patterns, all with sufficiently low background, and from structurally homogeneous samples.
Here, we present experimental results that address these problems and show the path towards single-particle imaging of macromolecules. We addressed the first challenges by aerosol injection of gold nanoparticles and achieved millions of patterns by using the megahertz-rate European XFEL [15] and a relatively large illumination area of the XFEL beam. The particles were chosen for the high scattering power of gold, which balances the reduced intensity from the large beam size to provide scattering signals at the levels expected from biological materials once tight focusing is achieved.
The challenge of structural heterogeneity is addressed computationally. Even though individual diffraction patterns contained as few as 0.0012 photons per pixel on average, we show that this is sufficient to not only extract the orientations of particles, but also to disentangle structural variations. We obtain a 3D structure approaching 2 nm resolution, which is significantly improved compared to what could be achieved without structural sorting.
With further improvements in aerosol sample delivery to increase the particle density in the X-ray focus [16][17][18], more highly focussed X-ray beams can be used to obtain similar data from biomolecules. The computational techniques developed here also open the way to experiments that can reveal thermodynamically rare states in an ensemble and characterise heterogeneous ensembles with statistical rigour. The short exposure times set by the femtosecond pulse duration will also offer unprecedented opportunities for capturing the dynamics of macromolecules in real time.   Figure 1: Experimental setup. XFEL pulses were focused by a series of Kirkpatrick-Baez mirrors into a 3×3 µm 2 spot and scattered off particles in the aerosol stream to produce diffraction patterns on the AGIPD. The lower inset shows the timing structure of the XFEL pulses at the instrument while the top inset shows representative SEM images of the cub42 and oct30 samples; scale bars are 100 nm. The low-resolution part of the detector used for the structural sorting is highlighted in green.

Diffraction data collection
Data was collected at the SPB/SFX (single particles, biomolecules and clusters/serial femtosecond crystallography) instrument [19] of the European XFEL using 6 keV photons focused into a 3×3 µm 2 spot, as measured by a 20 µm-thick YAG screen in the focal plane. Individual X-ray pulses were generated with 2.5 mJ of energy on average (2.6 × 10 12 photons). The pulses were delivered in 150-pulse trains with an intra-train repetition rate of 1.1 MHz and trains arriving every 0.1 s, leading to a maximum data collection rate of 1500 frames/second. A detector built specifically for this burst mode operation, the AGIPD [20], was placed 705 mm downstream of the interaction region to collect the diffraction patterns for each pulse individually up to a scattering angle of 8.3 • at the center-edge of the detector (see Fig. 1).
Gold octahedra and cubes, each of two different sizes, were sequentially injected into the X-ray beam using an electrospray-ionisation aerodynamic-lens-stack sample delivery system (see Methods). The nominal sizes of the particles measured using scanning electron microscopy were 30 and 40 nm for the octahedra and 42 and 17 nm for the cubes. In the rest of the article these samples are described using the codes oct30, oct40, cub42 and cub17 respectively. The octahedra and cubes were prepared using different protocols, generating different heterogeneity profiles as will be seen later.
Diffraction patterns were observed in around 10 % of the collected frames. This relatively high hit ratio compared to those achieved with biological particles in similar conditions was due to a combination of the relatively large X-ray focal spot size, high particle concentration and high mass and density of the larger gold nanoparticles, leading to lower speeds after acceleration by the gas flow in the aerodynamic lens stack [16,21,22]. Lower speeds lead to higher spatial densities, and thus higher hit ratios for the same particle beam size. Table 1 shows the statistics of the number of frames collected for each sample as well as the various filtration steps after the analyses described below.
When using the peak repetition rate of 1.1 MHz and 150 pulses per train, diffraction patterns corresponding to the shapes of cubes and octahedra could be observed, but a high fraction of the diffraction patterns appeared to originate from spherical particles (see Table 1 and third column of Fig. 2). This was found to be caused by the melting of particles in the wings of the previous XFEL pulse in the train, as the particles approached the focus. To reduce this occurrence we therefore reduced the intra-train repetition rate from 1.1 MHz to 550 kHz, providing only half the available pulses; further reduction of the repetition rate was tested but not found to be necessary. This reduced-rate mode was used to collect most of the data for the three larger samples (but not the cub17 sample).

Single hit selection by 2D classification
Frames with diffraction from particles were detected by setting a threshold on the number of pixels in the AGIPD detector that recorded at least one photon (see Methods). Unfortunately, not all the particles are of interest, even accounting for the heterogeneity. The extraneous patterns include those from spheres formed after melting, multi-particle aggregates and other possible contaminants. In previous work, either manual selection [10,13] or manifold learning methods [12,23,24] have been used to classify patterns and reject outliers. We adopt an alternative approach, similar to one commonly used in cryo-EM [25], but implemented in diffraction space. Two-dimensional orientation determination into multiple models was performed in the detector plane using the EMC algorithm [26,27] implemented in Dragonfly [28]. The in-plane rotation angle (θ) and relative incident fluence (φ) of each diffraction pattern was determined collectively and multiple independent 2D intensity models were reconstructed. Each of these intensities represent an average of aligned copies of a subset of the patterns from the whole set. In addition to the EMC algorithm being highly noise-tolerant [7,29,30], one can also use it to examine the average models to understand what type of particles are in the dataset.
In this experiment, 50 random white noise 2D intensity models were used as initial guesses to perform the classification for each sample, using only the low resolution part of the detector highlighted at this stage (see Fig. 1). Some of the reconstructed intensities are shown in Figure 2. The first two columns of the figure show representative examples of 'good' models of each sample, chosen manually to be those with high contrast and strong streaks for further processing. The third column shows an average of diffraction from rounded particles (except in the cub17 case where a dimer average is highlighted). These models were used to determine the sphere fraction shown in Table 1. Finally, the last column shows low-contrast models where a diverse set of particles were averaged.
The 2D classification also enabled the analysis of size-heterogeneity from those models where the faces of the nanoparticles were parallel to the X-ray beam. In these cases, one observes strong streaks on the detector and the fringe spacing indicates the distance between these parallel faces. The size distributions of the samples inferred this way are shown in Fig. 3(a). The octahedral samples had a much broader size distribution than the cubic ones. While some of the breadth of the peaks is due to apparent size variations when the faces are not being perfectly parallel to the beam, the much broader size distributions of the octahedra suggest that they had more heterogeneity.
In addition, the octahedra were also noticeably asymmetric, as seen in Figs. 3(c) and (d). These histograms were made by identifying patterns which belonged to models with two strong streaks (e.g. top left model in Fig. 2). Another run of 2D classification with just these two-streak patterns showed no variation in the angle between the streaks, but only in the fringe spacing. This is to be expected since the angle is fixed by the 111 growth    Due to the low polydispersity of the cubes, they were used to determine the incident fluence distribution of the X-ray beam. Since the Fourier transform of a cube is the product of three orthogonal sinc functions, the size fitting procedure also generated a predicted incident fluence. The distribution from 102 480 patterns is shown in Fig. 3(b), yielding a maximum fluence of around 60 µJ/µm 2 , which leads to a lower bound estimate of around 540 µJ in the focal spot from the measured spot size. The actual fluence was likely higher as the particles were not ideal cubes and the scattering efficiency is reduced at high fluences [31,32]. One can also see that most diffraction patterns were obtained with lower incident fluences, because the particles interacted with the outer regions of the X-ray focus.

3D reconstruction with structural sorting
The fraction of good hits used for 3D structure reconstruction varied from 17 % for the cube samples to around 60 % for the octahedra (see Table 1). The 3D intensity distribution was obtained using these patterns before recovering the structures by performing phase retrieval using the difference map algorithm [7,33]. For computational efficiency, the 3D orientations were first determined using the low-resolution part of the detector where the highest resolution was 3.3 nm. A refinement procedure similar to that developed for serial crystallography [34] was used with the whole detector to get the full-resolution 3D intensities. In this procedure, only orientations in the neighbourhood of the most likely orientation of a given pattern from the low-resolution run were searched.
The intensities recovered in this manner had noticeably lower contrast than the equivalent slices in the 2D models. From the size distributions seen in Fig. 3, this could be attributed to structural heterogeneity. To counter this, the patterns were probabilistically partitioned into five intensity volumes in a manner equivalent to the 2D classification procedure. However, the initial guesses were not random white noise, but rather isotropically stretched/scaled versions of the average models reconstructed above. Five models, with stretch factors ranging from 0.9 to 1.1 were used as these initial seeds. The rest of the reconstruction proceeded without any restraints between these models or any symmetry constraints.
Once again, this structural sorting was performed at low resolution before refining the orientations of a subset of patterns from a single model to get full-resolution intensities.
A comparison of orthogonal slices through the 3D intensity for the oct30 sample is shown in Fig. 4(a). The left column, showing the single model reconstruction with 1.4 million patterns has noticeably worse fringe contrast and background than the equivalent slices in the right column or in the first two columns of the 2D classification output shown in Fig. 2. The homogeneous set had 0.53 million patterns selected using the multi-model EMC reconstruction. The visual improvement is accompanied by an increase in the likelihood of the model intensities outside the central speckle for the common patterns in both sets, as shown in Fig. 4(b). The filled histogram shows the distribution of the per-pattern increase in likelihood, which we refer to as likelihood gain, while the two traces show the distributions for weak (relative scale 0.5 ± 0.1) and strong (relative scale 2.0 ± 0.1) patterns. The latter shows how brighter patterns are more selective towards an improved model. Figure 4(c) shows the same information for the oct40 sample, where the gain ratio is smaller, but still greater than 1. The 2D size distributions shown in Fig. 3(c) were re-calculated for each subset of patterns belonging to the five models and plotted in Fig. 4(d), confirming the different sizes for each model, but also exhibiting a simpler structure than that of the full dataset.
For the cubic particles, a single model 3D reconstruction was deemed sufficient, due to the relative monodispersity of the sample. The selection of 'good' hits from the 2D classification was more stringent, including only high-contrast cube-like patterns. The incident fluence factors were estimated in the first few iterations where the calculated probability distributions were broad and then later kept fixed (see Methods).
The electron densities were reconstructed by performing 3D iterative phase retrieval on the full-resolution intensity volumes (see Methods for details and Supplementary Fig. 1 for intensity slices). Figure 5(a) shows the reconstructed electron densities as isosurface plots. The contour levels were chosen where the gradient of the density was highest. The phase retrieval transfer function (PRTF) metric as a function of wavevector q is shown in Fig. 5(b). This metric is a measure of the reproducibility of recovered phases when starting from 128 random models. The 3D PRTF distribution was smoothed using a Gaussian kernel with a width equal to 1/3rd of the fringe width. The shaded region around each line shows the range of values in each q shell, highlighting the strong anisotropy of the metric due to the faceted nature of the objects. The intersection with the common 1/e threshold determining the resolution is shown in Table 1. The resolution normal to the flat faces is 2 nm or better for all samples, while the resolution is relatively low far from any strong streaks in Fourier space. This angle-dependent resolution is a property of the diffractive-domain averaging before phase retrieval, but also due to the strongly faceted shape and lack of internal structure of these objects, both of which are not representative of biological objects.

Discussion
We have demonstrated an order-of-magnitude increase in data collection efficiency along with much higher imaging resolution than previously achieved for X-ray single particle diffractive imaging, setting a template for future SPI experiments at the European XFEL and elsewhere. We have also shown that with these large data sets, one can structurally sort the particles and average a narrow size and shape range to obtain higher resolution. A similar problem is expected to be faced when imaging biological particles and the method developed here shows the way towards overcoming conformational variability in the Fourier domain.
Although we benefited from the strong scattering cross section of gold compared to organic materials, with the commissioning of a sub-micron focus at the SPB/SFX instrument, we can expect comparable signal strengths from organic materials. Unfortunately, smaller X-ray foci would also mean lower hit ratios with the current sample delivery setup. Improvements could be made through optimised focussing for the targeted size distribution [16] or cryogenic injections systems [18] which additionally allow conformational selection [35]. Another approach is to keep using the larger focus and conjugate the particles with gold nanoparticles to assist hitfinding and orientation determination [36]. The effective hit rate can also be increased by using more pulses from the European XFEL (max. 2700) than the AGIPD detector can save (max. 352) and vetoing in real time those frames which do not contain diffraction signal.
The class of experiments exemplified here can also be applied to study rare events such as transient states in a spontaneous phase transition or high free-energy states. Since each image is collected serially, one can identify relevant subsets corresponding to interesting states without averaging over all patterns. In this work, we have taken the approach of treating the objects as general 3D contrast functions with no a priori information.
One can also envision a parameterised refinement approach which should enable a finer characterisation of the structural landscape of the ensemble.

Sample preparation
The octahedral gold nanoparticles were synthesised using published protocols [37,38] with poly(diallyldimethylammonium) chloride (PDDA) polymer coating to avoid aggregation. The cubic particles were synthesized in water using the method described in Park et al. [39] with cetyltrimethylammonium chloride (CTAC) as a stabilising agent. In order to obtain the requisite 10 12 − 10 13 particles/mL concentration to approach an average of one particle per electrospray droplet [17], all syntheses were concentrated from initial values of 10 9 particles/mL for the cubes and 10 11 particles/mL for the octahedra and excess ligands were removed by centrifugation. Scanning and transmission electron microscopy images of the samples are shown in Supplementary Figs. 2 and 3.

Aerosol sample delivery
The samples were suspended in 10 mM ammonium acetate and aerosolized using an electrospray nebulizer (average flow rate 200 nL/min) and neutralized before delivering into the X-ray interaction point using the aerodynamic lens stack [17]. An electrospray differential mobility analysis (ES-DMA) setup was installed at the beamline and the particles generated by the electrospray could be diverted into the ES-DMA to characterise the size distribution and concentration of aerosolized particles. Particle size distribution measurements were carried out with an electrostatic classifier (TSI 3082) together with the DMA (TSI 3081). The DMA was connected to a condensation particle counter (CPC, TSI 3789). Representative size distributions are shown in Supplementary Fig. 4. To diagnose the width and density of the particle stream in the X-ray interaction region after aerodynamic focusing, we employed a Rayleigh scattering diagnostic (see Supplementary  Fig. 5) using a frequency-doubled Nd:YAG laser which was mirror-incoupled perpendicular to both the X-ray beam and particle stream [21,22]. These two diagnostic tools helped in assessing both the quality of the samples as well as the transmission efficiency of the sample delivery system.

Online monitoring
In order to help align the experiment and dianose problems during data collection, the Hummingbird software [40] was connected to the Karabo bridge in the European XFEL DAQ system [41] to receive data with a delay of a few seconds. Since most of the photons in an SPI experiment are concentrated at low resolution, only the module of the AGIPD closest to the beam centre was used for online analysis. The use of uncalibrated data and only a single module enabled a frame rate of up to 800 Hz using all 176 memory cells of each pixel of the AGIPD available in this experiment. The analyses conducted live included lit-pixel hit finding and hit-ratio determination (see the following Preliminary analysis section), sphere model size determination and the detection of the fraction of spherical particles by analysing the azimuthal variation in intensities. The latter was used to understand and fix the particle melting issue mentioned in the Results section.

Preliminary analysis
The AGIPD detector was calibrated using offset constants for each cell in each pixel.
Except for a few pixels near the beam centre, no pixels switched gain mode. After offset correction, the number of pixels containing at least 0.7 of a photon was calculated for each frame. In the absence of particles, the number of such pixels is normally distributed with a mean dependent on background from the beamline, carrier gas and detector false positives. A threshold of 3σ over the mean number in each run was used to select frames with particle scattering. Over the entire experiment, the hit ratio fluctuated between 7%-15%. For this and future analyses, memory cells in pixels with outlier dark offsets or dark noise were masked out, along with the double-wide pixels along ASIC edges.
These hits were converted to photons by first subtracting the dark offsets, correcting for per frame common mode shifts by subtracting the median of each 64x64 pixel ASIC and then subtracting a pixel-wise running median of the last 128 frames over all cells. The last step was important in removing artifacts due to the slow drift of dark offsets on a pixel level. These corrected detector values were then converted to integer photon counts by thresholding with a variable cutoff using the following procedure.
The probability distribution of detector ADUs (analog-to-digital units) at a pixel in the absence of photons is a Gaussian centered at 0 with a cell-dependent width. The 1-photon distribution is a shifted copy of the 0-photon distribution with a height which depends on the signal level (ignoring charge sharing for the large 200 µm pixels). The optimal threshold was chosen to be the point at which these 0-photon and 1-photon distributions intersect for a signal level of 10 −3 photons/pixel. If the 1-photon distribution is centered at m 1 ADUs and the standard deviation of the noise of a cell is σ ADUs, the threshold This threshold minimises the total error rate (false positive plus false negative) due to detector noise at the chosen signal level, ignoring charge sharing effects which are small for the 200 µm pixels of the AGIPD. For higher signals, at lower resolution, the error rate would be higher, but biased towards false negatives. See Supplementary Fig.6 for the effects of the various corrections on the integrated detector image. For the current detector configuration and photon energy the 1-photon peak was centered at 47 ADUs and the average threshold was 0.755 of a photon.

Intensity reconstruction
Both two-and three-dimensional intensity volumes were reconstructed using the Dragonfly package. Detector files were generated and refined manually starting from initial geometries from a previous serial crystallography experiment [42]. The geometry refinement only involved adjusting the positions of the detector quadrants since the modules within a quadrant had not been moved between the experiments. Two detector files were produced, one for the inner eight 128x64 pixel detector ASICs, while a high-resolution version contained all 1024x1024 pixels.
The photon-converted data was saved in the sparse Dragonfly .emc format with file sizes of around 10 GB per sample. For the 3D reconstructions, the detector files specify the reciprocal space voxel coordinate of each pixel, which involves defining the radius of curvature of the Ewald sphere in voxels. The natural choice of the detector distance in pixel units (3525) produced too large an oversampling factor, with a fringe spacing of around 20 voxels for the largest cub42 sample and even higher for the smaller samples. For the low-resolution detector file, the radius of curvature was set to be 2000 voxels, generating 253 3 voxel volumes. For the full-resolution detector, it was 1500 voxels for all samples except cub17, where it was 1000 voxels.
For all reconstructions, the deterministic annealing procedure [28,43] was used to improve the convergence of the algorithm. The annealing parameter β d was initially set for each pattern, d based on the number of scattered photons using the following empirical formula: where C d refers to the number of orientationally relevant photons in the pattern. This generates a lower value for brighter patterns, broadening their otherwise sharp probability distribution over orientations. The specific expression was tested in simulations to produce a relatively flat dependence of the mutual information I(K, Ω) [28] on the signal strength. The parameter was increased by a factor of 2 every 10 iterations for each pattern.
In the 2D reconstructions, 180 angular samples were chosen in the range from 0 to 2π for each model. The low resolution 3D reconstructions were performed with an orientational sampling level of 8, which corresponds to 25 680 samples [26]. The high-resolution refinement went up to a sampling level of 20 (400 200 samples). For the 3D multi-model reconstructions, the initial intensities were generated by isotropically stretching the singlemode intensity volume using linear interpolation and the initial fluence factors were also used from that run.
For the cubic samples, the 3D reconstruction pipeline was modified. First, 40 iterations were performed with the initial β d parameters without an annealing schedule. For the larger cub42 sample, only the pixels corresponding to the first 3 diffraction fringes (12.9 nm resolution) were used to determine the orientations and fluence factors. This was to avoid instabilities since most of the Fourier power at higher resolution is concentrated in the streaks normal to the faces and the angle between the streaks is large enough that interference between them poorly constrains the model. Once the low resolution, rotationally blurred intensities were stable, the fluence factors were fixed and the annealing schedule was enabled. The rest of the reconstruction proceeded in a similar manner to the octahedra. Isosurface plots for the three larger particles using the low-resolution data are shown in Supplementary Movie 1.

Size and incident fluence fitting
The 2D classification was first used to identify patterns where strong streaks were visible. These classes can be seen in the first column in Fig. 2. A pair of parallel faces on a particle produce a sinc-function dependence in the Fourier transform along the face normal. This was used to determine the interfacial distance from each pattern.
The classification procedure not only allowed us to identify the patterns which have strong streaks, but was also used to determine the angles of these streaks since we knew by what angle the pattern had to be rotated to fit the model. For these selected patterns, the intensity distribution along the streak was calculated by integrating over a 21-pixel wide strip along the streak. The size was determined by cross-correlating the intensity distributions with sinc functions for sizes from 10 to 50 nm in 0.1 nm increments. The size was chosen as the one which produced the maximum Pearson correlation coefficient and only those streaks with a coefficient greater than 0.9 were included. Figures 3(c, d) were generated by only considering patterns with 2 strong streaks.
For the cub42 sample, the brightness of the two-streak patterns were used to determine the incident fluence by assuming that they were generated by perfect cubes. A procedure similar to the sphere-sizing done in previous works [44] was used to determine the incident fluence using the scattering cross-section for gold at 6 keV.

Phase retrieval
Before performing phase retrieval, the intensity volumes were processed in the following manner: Background subtraction was performed using a rolling minimum filter with the window size 1.5 times the width of a fringe using the ndimage.minimum filter in SciPy [45]. Since no fringe contrast was visible at the outer resolution edges (see Supplementary Figure 1), the data was truncated such that the corners of the cube were within the sphere. The full-period resolution of the cube at the center-edge was 1.61 nm and there were 384 3 voxels for the three larger samples and 256 3 voxels for the cub17 sample.
A combination of the error reduction (ER) algorithm and the difference map (DM) algorithm [33] were used to reconstruct the electron densities from the background-subtracted intensity distribution. For each phasing run, 400 iterations were performed, 100 ER, 200 DM and 100 ER. Within each iteration, a dynamic real-space support constraint was applied by sorting the electron density values and only keeping the top N supp values. The volume N supp was chosen such that the histogram of densities inside the support had a small fraction of low values to ensure that the support was not too tight.
128 phasing runs from random white noise initial guesses were performed for each sample and the resulting densities were aligned and averaged to produce the final electron densities as well as the phase retrieval transfer function (PRTF). The PRTF strongly depends on the intensity at a voxel and thus exhibits an oscillatory behaviour along the fringes, which is not reflective of the quality of the structure at a given resolution since the lack of accurate phases near an interference minima barely affect the real-space structure. Thus, the 3D PRTF distribution was smoothed with a Gaussian kernel with a width half that of a fringe [13]. The 3D isosurface plots were rendered using Chimera [46] and the contour levels were determined by using the "Surface Color" feature, colouring the surface by the gradient of the density and choosing the level which had the highest density gradients.
with contributions from all authors.  The buffer peaks near 10 nm are absent after passing through the aerodynamic lens to reach the interaction region and will, in any case, not be focused in the same region as the heavier AuNPs. Note that the ES-DMA measures the electrical mobility diameter which is dependent on both the size and shape of the particle [47]. Top: Histogram of detected particle positions in the interaction region for the cub42 sample. The histogram was obtained from 1000 trains with one image per train. The X-ray beam propagates along the Z-axis and the Y-axis is vertical in the laboratory frame. The finite length of the particle beam in the vertical direction is due to the ∼250 µm optical laser spot size. The particle positions were detected using the peak local max function of Scikit Image [48] after appropriate background correction. Bottom row : Three representative single backgroundcorrected images from the same run.  Figure 6: Detector corrections. Integrated photon-converted detector patterns for all 44 800 hits in a single run with progressively more detector corrections. First only standard dark offset and common mode corrections were applied before thresholding at 0.7 of a photon. Next, the running median offset value over the last 128 trains was subtracted from each pixel before thresholding. Finally, a variable threshold was used for each memory cell depending upon the standard deviation (sigma) of the cell. All images have the same colour scale, saturating at 10 −3 photons/pixel/frame. The dark gray background shows panel gaps and masked pixels. The fourth plot shows the radial average intensity for the three images.  Fig. 5(a).