HOTCAP: a new software package for high-speed oscilloscope-based three-dimensional bunch charge and position measurement

A software package to be used in high-speed oscilloscope-based three-dimensional bunch-by-bunch charge and position measurement is presented. The software package takes the pick-up electrode signal waveform recorded by the high-speed oscilloscope as input, and it calculates and outputs the bunch-by-bunch charge and position. In addition to enabling a three-dimensional observation of the motion of each passing bunch on all beam position monitor pick-up electrodes, it offers many additional features such as injection analysis, bunch response function reconstruction, and turn-by-turn beam analysis. The software package has an easy-to-understand graphical user interface and convenient interactive operation, which has been verified on the Windows 10 system.


Introduction
New-generation synchrotron radiation light sources have the following basic characteristics: ultra-low beam emittance (a few nmrad or even tens of pmrad), ultra-high beam stability (with orbital feedback control accuracy being mostly at the micrometer or submicrometer level), high average beam current ( [ 300 mA), a small-aperture beam vacuum pipe (\30 mm in diameter), a large number of vacuum inserts, top-up operation mode, small dynamic aperture, and high time resolution experimental ability (of the order of picoseconds) [1,2]. The above-mentioned characteristics make how to minimize instability a key technological problem to solve. New-generation light sources have higher requirements for beam stability. Accurate beam diagnosis technology guarantees smooth operation of synchrotron radiation facilities. Bunch-bybunch measurement technology is a popular research direction in the field of beam diagnosis. Using bunch-bybunch measurement technology enables accurate beam diagnosis for each bunch rather than just for the average status of all bunches. This is beneficial for the operating staff to better understand the operating status of the device, facilitating targeted optimization to ensure that the synchrotron radiation device has a better light supply quality [3][4][5].
Developing a set of bunch-by-bunch signal processing algorithms and establishing a bunch-by-bunch parameter measurement system with high resolution and high refresh rate are crucial. In the early days, limited by electronics technology, beam measurement generally performed on the storage ring mainly included beam closed-track measurement [6] and turn-by-turn beam measurement [7,8]. In recent years, with the development of electronic technology and the need for accelerator machine research, high-precision bunch-by-bunch measurement has become important in accelerator research around the world [9][10][11].
Bunch-by-bunch measurement mainly includes transverse position measurements [12][13][14], longitudinal phase measurements [15][16][17], and charge measurements [18][19][20], as well as beam size measurement and emittance [21][22][23]. However, position measurement is the most important and widely used method. The commonly used measurement methods mainly include two categories: phase-locked sampling with an external clock and random phase sampling with a constant frequency. The phase-locked sampling method with an external clock is usually based on a commercial bunch-by-bunch feedback system or a highspeed acquisition board. The random phase sampling method with a constant frequency is usually based on a digital oscilloscope.
Compared with high-speed acquisition boards and commercial feedback systems, digital oscilloscopes usually have a higher bandwidth and sampling rate. The random phase sampling method can be more convenient for performing secondary development to further study the bunches at different accelerators [24][25][26]. This makes it possible to build a bunch-by-bunch system that is common to all major accelerators worldwide. Therefore, it is advantageous to use a digital oscilloscope as a signal acquisition device for bunch-by-bunch measurement [27][28][29].
The sampling mode of a digital oscilloscope is usually asynchronous; that is, the signal acquisition does not use the accelerator's radio-frequency (RF) signal as an external clock. Signal acquisition is completely triggered by an internal clock. The advantage of this is that all acquisition results and further processing results are based on the absolute time when the accuracy of the internal clock is much better than the measurement accuracy. However, the high sampling rate and the non-phase-locked sampling mode result in a large amount of input data and complex digital signal processing. Currently, there are no comprehensive data processing algorithms or corresponding data processing software packages that are widely used in highspeed oscilloscope-based three-dimensional (3D) bunch charge and position measurements.
Based on the cross-correlation method, the Shanghai Synchrotron Radiation Facility (SSRF) Beam Instrument Group has studied a set of signal processing algorithms for high-speed oscilloscope-based 3D bunch charge and position measurement and released a visualization software package called HOTCAP. Each step of the parameter setting, data processing, and result reporting can be fulfilled using a mouse click.
The theoretical basis for this package is briefly described in Sect. 2. The main features and capabilities of HOTCAP are presented in Sect. 3. To further verify its usefulness, test experiments at the SSRF and at the Hefei Light Source (HLS) are presented in Sect. 4. Section 5 provides a summary of the results.

Theoretical basis
The button beam position monitor (BPM) is usually used to measure the position of the beam. There were four button-type electrodes within the BPM. The signal obtained from the button-type electrode contains almost all the bunch-by-bunch information.
When a bunch passes through the button electrode, the electrical signal is drawn from the electrode by the feedthrough, which contains almost all the parameter information of the bunch. The electrical signal collected by the collection device can be expressed as where I(t) is the current intensity of the resistor. A bunch with a longitudinal Gaussian distribution can be expressed in the time domain as I 0 ðtÞ [30,31]. It can be seen that the collected voltage signal is related to the transverse position (Fðd; hÞ), longitudinal position (longitudinal phase t 0 ), and charge (e Á N) of the bunch. Unfortunately, the actual collected signal is not completely in accordance with the theoretical derivation, and the limited bandwidth of the feedthrough introduces time delays. In addition, the bandwidth limitation of the acquisition device also changes the signal waveform in the time domain.
Therefore, it is theoretically feasible to build a bunchby-bunch 3D charge system by using the button BPM, which is common in almost all accelerator storage rings.

Main features and capabilities
In this system, button BPMs are used as the bunch signal pick-up. A high-speed broadband oscilloscope was chosen to record the signals from the four button electrodes of the BPMs. These signals can be directly read into the memory using the software package HOTCAP as input to HOTCAP and processed to obtain various physical parameters that characterize the beam current and machine state. The system framework is shown in Fig. 1.
The software package was developed based on Python, and the visual interface was based on PyQT5. To improve the operation speed, NUMBA was used as the acceleration scheme of the compiled script in functions with a large amount of calculations. The software package was divided into three parts: the user interface (UI) module, the calculation module, and input-output (IO) module. The calculation module was run in a separate child thread.
The software package performs bunch-by-bunch 3D parameter extraction, bunch charge extraction, refilled charge parameter extraction, and analysis of the BPM local wake field. The parameter calculation function for this part was realized in the calculation module. The IO module acts as the data interactive interface between the software package and MATLAB, which reads the original oscilloscope file and saves the calculation results into the corresponding format output file. The UI module serves as the console for the user. Although it reads the machine parameters provided by the user as the necessary parameters for calculation, it shows the calculation results to the user in different forms (e.g., bunch-by-bunch-by-turn, filling analysis, etc.). A simplified schema of the system is presented in Fig. 2.

UI module
The user's interaction with HOTCAP is entirely based on an interactive front-end interface. The interface was built with Qt Designer, and its underlying framework was developed based on PyQT5. At present, all user operations can be completed by clicking the mouse. The user must select the data file to be calculated in the console and select the calculation method. By inputting the necessary machine parameters and calculation mode through the interface, HOTCAP sends the data to the calculation module for the corresponding calculation. After the calculation is completed, the calculation results can be displayed in different forms in the UI (bunch-by-bunch display, turn-by-turn display, and refilled charge movement pattern display).

IO module
The IO module realizes the storage of calculation results and data interaction between the software package and the workstation. Currently, HOTCAP's input and output file format is mat. The mat file can be viewed using MATLAB. Such files usually exist in MATLAB MAT-file format. The MAT file extensions are primarily grouped under the data file category. In less common applications, they can also be uncommon files or 3D image files. MAT files have already been supported on Windows, Mac, and Linux platforms. They are compatible with desktop computers and mobile computers. HOTCAP is compatible with Àv7 and lower versions of mat-type files.
The package accepts two or four channels of oscilloscope data, which can be collected from two pairs of electrodes in one direction or from four electrodes in two directions. The current IO module requires these data to be stored as a mat file with variable names of BPM1-BPM4. The calculated results can be stored in the workstation in mat form; these results include but are not limited to the 3D position and charge of all bunches, the response function, the current RF, and the local wake field.
In addition, the IO module also undertakes the function of storing the calculation results. When the calculation module has completed the calculation of one set of data, the results are stored in the IO module, freeing the calculation module to calculate the next set of data. The stored calculation results can be displayed on the interface in the form of line charts, bar charts, etc.

Calculation module
The calculation module is the core module of the HOTCAP software package. It is responsible for almost all signal preprocessing and bunch parameter extraction functions. The calculation module is composed of a main calculation function and 17 subfunctions.
The following section introduces the calculation process and the core algorithm of the calculation module in detail. The calculation process is shown in Fig. 3. The primary core algorithms are used for longitudinal phase extraction based on cross-correlation, response function reconstruction based on software resampling, and refilled bunch signal extraction [32].

Preprocess
Data preprocessing mainly includes normalization, DC compensation, signal initial peak finding, and RF coarse and fine adjustments.  Normalization eliminates the influence of different sampling equipment and normalizes the amplitude of the signal to a fixed range. The DC compensation function removes the DC bias caused by the acquisition equipment by recording long-term statistics of the value at the null signal (empty bucket). The function of the signal initial peak finding function is to find the peak value of the first bunch at the head of the first bunch train, which is used as the starting point for signal processing.
RF adjustment is the most important part of data preprocessing. Although the RF for a storage ring is usually fixed, in reality, it constantly changes slightly. This RF shift, though small, is very serious for precise 3D bunchby-bunch measurements. Therefore, it is necessary to adjust the RF in real time according to the bunch signal.
The RF adjustment includes two steps: coarse and fine adjustments. During coarse adjustment, the peak point of the signal collected during each turn of the arbitrarily selected probe bunch is used as the feature point, and the time interval between the feature points is calculated. Using the time interval, the RF is then adjusted coarsely. The fine adjustment of the RF is more complicated, and the longitudinal phase of the probe bunch needs to be calculated first. The standard of the fine adjustment is that the equilibrium value of the longitudinal phase remains unchanged. The calculation method for the longitudinal phase is introduced in the following subsections.

Response function reconstruction
The response function described in this section refers to the shape of the signal pulse visible to the acquisition device when the bunch passes through the electrode of the BPM, that is, the response of the electrode and the acquisition device to the bunch.
The reconstruction of the response function recovers from the signal the shape of the signal of each bunch, that is, V(t), with a particular bandwidth. The signal excited by a bunch passing through the BPM collected by a broadband device is very short-lived ( $ 1 ns at the SSRF), which means that the acquisition device can only collect a few sampling points. The signal actually has several points on the continuous smooth response function. However, what is needed is a continuous and smooth response function with a higher sampling rate. The reconstruction response function uses multiturn data to splice a response function with a high sampling rate.
For almost all storage rings, the sampling frequency f of the oscilloscope must adhere to the following formula: where m is a positive integer and f c is the cyclotron revolution frequency. When the bunch passes through the BPM again, the acquired signal point will not be in the same position as the response function last time. t 0 can be obtained when the bunch passes the BPM every turn because f c and f are known. Based on removing the initial phase t 0 , through the splicing of multiturn sampling data, a high sampling rate response function can be constructed.
Owing to the signal noise and oscillation of the transverse position of the bunch, the high sampling rate response function constructed by multiturn splicing is not very smooth. From the perspective of the frequency domain, the splicing signal is mixed with high-frequency noise. From the perspective of the time domain, the line of the spliced response function will be thick. The real response function is the most probable position of the response function obtained by direct splicing. Therefore, this computing module achieves the effect of low-pass filtering through window smoothing to obtain the best response function. A random example was selected to show the process and effect of the reconstruction (shown in Fig. 4).
In HOTCAP's built-in module, the equivalent sampling rate of the reconstructed response function is 100 THz; that is, the interval between adjacent points is 0.1 ps. If a higher measurement accuracy (of \0:1 ps) is required, the equivalent sampling rate will need to be increased. This part allows parameters to be adjusted in the secondary development of the HOTCAP package. Because the influence of the measurement error of common sampling devices is often [ 0:1 ps, the default value of the equivalent sampling rate is set to 100 THz, which helps to reduce the amount of computation.
Because the response functions of different bunches collected by different electrodes are different, they are all reconstructed independently. In addition, because the response function is obtained by the splicing of multiple measurement results, this algorithm acts on the premise that the response function of the bunch does not change dramatically in a short time. This requires that the bunch length and bunch charge do not change dramatically in a short time ( $ 1 ms). Most stable synchrotron radiation devices meet this requirement. If this condition is not met in a special machine status, the shape of the response function with different bunch lengths needs to be constructed and calibrated in advance. The HOTCAP package does not currently support this extremely special machine status, requiring secondary development of the package by the user.

Phase and position extracting method
The longitudinal phase and transverse amplitude extraction are the most important parameters in the HOT-CAP package. The longitudinal phase extraction is based on the sliding cross-correlation method. The transverse position is obtained by the difference ratio sum of the signal amplitudes of the four electrodes. The amplitude of the signal is then obtained by a linear fitting.
As shown in Fig. 5, the blue line represents the response function of the signal when a bunch passes through an electrode. The response function represented by the blue line is composed of a sampling point in thousands of turns. The yellow line represents the temporary response function excited when the bunch passes through the electrode in the Nth turn. The relationship between the blue standard response function and the yellow temporary response function in the Nth turn can be expressed as where Func is the response function with a blue line, amp N is the relative amplitude of the Nth turn, and u N is the relative longitudinal phase in the Nth turn. These are very important for the calculation of the transverse and longitudinal phases. In fact, the temporary response function represented by the yellow line is not available because the sampling rate of the acquisition devices is insufficient. We obtain a number of discrete sample points (represented by purple ? signs). The standard response function represented by the blue line is obtained with the response function reconstruction method presented in the previous section. Therefore, the main function of the algorithm module described in this section is to calculate the phase and amplitude with the response function represented by the blue line and the data points acquired in each turn represented by the purple ? signs as input data.
The data points acquired in each turn are actually a continuous and equally spaced sampling of the temporary response function in this turn, which can be expressed as where d n N is the data point array in the Nth turn, n is the sampling index, and t 0 is the initial sampling phase, which can be obtained using f and f c . Using the sliding crosscorrelation method, the relative phase is obtained by the sampling array fd n N g and the standard response function (Func). The operation method is described in detail in the following section.
The first step is to resample the standard response function using the oscilloscope sampling interval as the sampling interval. The sliding in cross-correlation means that software resampling is repeated several times in succession with different initial sampling phases (uðiÞ). The resampling can be expressed as where LUT n;i is the resampling data point array with the initial sampling phase (t i ). From Eqs.
With the Pearson correlation coefficient chosen as the correlation coefficient, the cross-correlation coefficient can be expressed as Mathematically, this coefficient represents the cosine of the angle between two vectors, ranging from -1 to 1. If and only if Eq. (6), the correlation coefficient is equal to one (amp N [ 0). In this case, u N can be obtained from Simultaneously, the relative amplitude (amp N ) in the Nth turn can be obtained by a linear fitting: Thus far, the phase (u N ) and relative amplitude (amp N ) in Eq. (3) are calculated. Clearly, u N is the longitudinal phase. The absolute amplitude with different electrodes (A, B, C, or D) is the relative amplitude multiplied by the amplitude of the corresponding standard function. The peak-to-peak value of the standard response function is taken as the amplitude of the standard response function. The bunch-by-bunch transverse position (x, y) is expressed as follows: where k x and k y are constants. Eq. (10) is valid only when the beam is at the center of the vacuum chamber. The offset beam includes a systematic error. For most synchrotron radiation facilities, the beam is very close to the center of the vacuum chamber, so this algorithm can meet the accuracy requirements. It should be implemented to include a polynomial fit to reproduce the monitor chamber structure if the beam is far from the center of the vacuum chamber [33,34]. This is not the core part of the software.
If the user has a special machine condition, a nonlinear correction can be added by modifying the HOTCAP code.

Refilled charge signal extracting method
The top-up operation mode requires frequent injection to maintain a constant high current. This top-up injection process is a perturbation to the storage ring, which deflects the beam far away from the equilibrium point. The refilled bunch contains stored and refilled charges. The signal collected by the BPMs is a combination of stored and refilled charges. During the injection process, the 3D movement of the newly refilled charge carries different information from the 3D movement of the stored charge [24].
The HOTCAP software package contains an algorithm module that extracts the refilled charge signal from the refilled bunch signal, enabling the evolution of the 3D position of the refilled charge to be further extracted. The principles of this module are briefly described in this section.
The refilled bunch consists of the stored charge and refilled charge that has just been injected. For a linear measurement system, the acquired signal of the refilled bunch can be expressed as signalðnÞ ¼ signal s ðnÞ þ signal r ðnÞ; ð11Þ which is the sum of the refilled charge signal (signal r ðnÞ) and the stored charge signal (signal s ðnÞ). To obtain the refilled charge signal, the stored charge signal needs to be subtracted from the acquired signal. However, it is very difficult to obtain the signal of the stored charge because it is mixed in the sum signal. The algorithm used in the HOTCAP software package is implemented as follows: The first step entails finding the injection transient and refilled bunch index by the sudden change in the charge. In the second step, the 3D position of the stored charge in the refilled bunch is obtained indirectly using the common mode of the 3D position of all bunches. The third step involves reconstructing the stored charge signal according to the 3D position of the stored charge. In the last step, the reconstructed stored charge signal is subtracted from the acquired signal to obtain the refilled charge signal.
In most synchrotron radiation facilities, adjacent bunches have similar 3D positions. Therefore, the average positions of several stored bunches (without refilled charge) near the refilled bunch were selected as the stored charge positions in the refilled bunch. For HOTCAP, there are 10 nearby bunches. The specific value should be selected according to the situation of different machines. Because it is necessary to know the charge of the refilled charge and the charge of the stored charge when rebuilding the signal, the response function of the stored charge is constructed in advance before the injection occurs.

Experiments at the SSRF
It is worth mentioning that, because the 3D position of bunches is constantly changing, it is necessary to calculate the 3D position of stored bunches at each moment (each turn) in advance. The purpose is to obtain the 3D position of the stored charge in the refilled bunch at each turn. This is very important for refilled charge signal extraction.

Test examples
To test the reliability of the developed package HOT-CAP, experiments were conducted at two different synchrotron radiation facilities: the SSRF and the HLS.
The SSRF, which consists of a 150-MeV linear accelerator, a 180-m, 3.5-GeV booster, and a 432-m, 3.5-GeV storage ring, is a multibunch, high-energy third-generation light source. The harmonic number was 720 at an RF of 499.654 MHz. The parameters of the SSRF are listed in Table 1. To improve the efficiency and quality of the light, a top-up filling pattern was adopted at the end of 2012, which resulted in more frequent beam injections and storage of $ 500 bunches with a 2-ns spacing [10,35].
Two different high-speed oscilloscopes, one with a sampling rate of 10 GHz and a bandwidth of 4.2 GHz and another with sampling rate of 25 GHz and a bandwidth of 10 GHz, were used, respectively, in the experiment. The 3D position evolution of the stationary light supply time and injection transient was studied.
Typical results plotted in this section were processed using the data recorded on July 25, 2020, unless otherwise noted (shown in Fig 6).
The results of typical injection transient record processing are shown in the package interface. In the experiment, the trigger of the oscilloscope was set to capture the injection transient, and the data of four channels (from four electrodes) were imported and loaded into the software package. After the calculation was completed, the results were displayed in four subinterfaces (parameter setting, bunch-by-bunch, turn-by-turn, and injection analysis).
As can be seen from the left figure of the parameter setting subinterface (Fig. 6a), the filling pattern of the SSRF is four trains, each consisting of more than 100 bunches with roughly equal charges. The 3D position, response function, and correlation coefficient of any bunch are shown in the bunch-by-bunch subinterface (Fig. 6b). As can be seen from the horizontal position shown on the upper left, the refilling process occurs around the 800th turn. Subsequently, the stored bunch is pushed far away from the equilibrium point by the impact of the kicker's field and displays a typical damped oscillation. By setting the bunch number, the user can display the information about the bunch that the user wants to study. In the turn-byturn subinterface (Fig. 6c), the weighted average of each bunch is weighted by the charge to obtain the 3D turn-byturn information. The instability caused by the filling process can be clearly seen in the turn-by-turn 3D position. The 3D position of the refilled charge is shown in the injection analysis subinterface (Fig. 6d). After the completion of the refilled charge extraction, the refilled charge with the maximum charge was selected as the calculation object. Its 3D position parameters and charge were calculated and displayed in the subinterface. From the longitudinal phase figure on the upper right, it can be clearly seen that, after injection, the longitudinal phase of the refilled charge exhibits a typical longitudinal damping oscillation, which is consistent with the theory.
The HLS, which consists of a linear accelerator and a 0.8-GeV storage ring, is a multibunch, high-energy light source. The harmonic number is 45 at an RF of 204.03 MHz [36]. The parameters of the HLS are listed in Table 1.
To improve the efficiency and quality of the light, a top-up filling pattern was adopted , which resulted in more frequent beam injections and storage of $ 35 bunches with a spacing of 5 ns. It can be seen that the machine parameters of the HLS and SSRF are quite different, making them suitable choices to verify the good generalization of HOTCAP.
The HLS is a second-generation synchrotron radiation accelerator built in the last century. Compared with the third-generation synchrotron radiation accelerator, the beam 3D position is not stable and the difference between the bunches is very large; that is, the 3D position of the adjacent bunches is very different. We analyzed the possible reasons for this phenomenon. First, there is a roomtemperature accelerating cavity at the HLS but not a superconducting accelerating cavity. In addition, highorder harmonic cavities may also increase the difference between bunches. We hope to use the HOTCAP software package to study the instability phenomenon and the difference between bunches at the HLS.
HOTCAP supports exporting bunch-by-bunch results in a file. Three typical bunches (a, b, and c) were selected from the 35 bunches stored in the HLS storage ring, and their 3D position and Pearson correlation coefficient evolutionary process are shown in Fig 7. It can be seen that their 3D positions at the same time are different, which proves the feasibility of the bunch-by-bunch measurement; that is, different bunch statuses can be measured simultaneously. In addition, the bunches in the HLS have a longitudinal phase oscillation. The oscillation amplitude was $ 100 ps, and the oscillation amplitudes of different numbered bunches were different at the same time. We suspect that this phenomenon is caused by the room-temperature accelerating cavity and the high-order harmonic cavity.
It is interesting to note that the amplitude of the horizontal position (red lines in a, b, and c) and the longitudinal phase oscillations (black lines in a, b, and c) are consistent,  It can be seen from the blue line (bottom figure) of each bunch that the correlation coefficient is no longer stable near 1 but is often \1, which indicates that the effect of cross-correlation matching is not perfect. The reason for this problem is that the bunch length of the HLS is not constant but changes greatly in a short time, which leads to a constant change in the response function. HOTCAP uses a fixed response function as matching data, which affects the calculation accuracy. The correlation coefficient decreases (in $ 30,000 turns) when severe instability occurs. This indicates that the bunch length in the HLS storage ring becomes long and short during this time.

Conclusion
3D bunch-by-bunch position measurements are very important for beam diagnosis. At present, there is no mature online measurement system that can simultaneously measure the bunch-by-bunch 3D position. In addition, the injection transient is an important window for machine research. Many synchrotron radiation devices operate in the top-up mode, which means frequent refilling. At present, there is no mature online system that can strip the refilled charge signal. Therefore, the Beam Instrument Group of the SSRF developed a 3D position measurement package called HOTCAP. HOTCAP can measure the 3D position of each bunch synchronously and strip the refilled charge signal to further analyze the injection transient. In this study, using a signal acquisition device with a digital oscilloscope, we realized an online measurement system that can simultaneously measure the bunch-by-bunch 3D position and charge. Depending on the developed versatile and user-friendly software package, it also offers many advanced features such as injection transient analysis, bunch response function reconstruction, and turn-by-turn beam analysis. This work is expected to be widely used in similar accelerators and can be further improved.
The HOTCAP package is highly versatile and can be applied to most synchrotron radiation facilities. At the same time, the current version of the software package requires that the bunch length in the storage ring be roughly constant. The reason for this requirement is that, if the bunch length changes dramatically, the response function will change accordingly, and the accuracy of the phasematching method based on cross-correlation will be reduced. Given that the bunch length is constant in most third-generation synchrotron radiation facilities, HOTCAP can meet most usage requirements. This problem can be solved in the future by upgrading the algorithm. The specific method would be to scan and calibrate the response functions under different bunch lengths for algorithm matching after actively changing the set value of the bunch length. With the help of this algorithm, HOTCAP would no longer require a constant bunch length but could also calculate the current bunch length in real time from the matching situation.
The current version of HOTCAP is available on GITHUB. 1 Peers with relevant research needs are free to download and use the package. In the future, we plan to open the source code to the platform so that other peers can redevelop it for specific accelerators.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.
Author Contributions All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Xing-Yi Xu, Yong-Bin Leng, and Yi-Mei Zhou.