Cosmic Ray Removal in Fiber Spectroscopic Image

TL;DR

This paper introduces a novel method for detecting and removing cosmic rays in single-exposure, multi-fiber spectroscopic images, achieving high detection accuracy and effective correction.

Contribution

The paper presents a new cosmic ray removal technique specifically designed for single-exposure, multi-fiber spectra using 2D profile fitting and a position-dependent noise model.

Findings

Detects up to 80% of cosmic rays

Corrects cosmic ray pixels with 97.8% accuracy

Works effectively on simulated and real data

Abstract

Single-exposure spectra in large spectral surveys are valuable for time domain studies such as stellar variability, but there is no available method to eliminate cosmic rays for single-exposure, multi-fiber spectral images. In this paper, we describe a new method to detect and remove cosmic rays in multi-fiber spectroscopic single exposures. Through the use of two-dimensional profile fitting and a noise model that considers the position-dependent errors, we successfully detect as many as 80% of the cosmic rays and correct the cosmic ray polluted pixels to an average accuracy of 97.8%. Multiple tests and comparisons with both simulated data and real LAMOST data show that the method works properly in detection rate, false detection rate, and validity of cosmic ray correction.

Tables

Table 1. Table 1: Results of CR Detection on Simulated Images

Items	IMG600	IMG1800
CR added	227451	227451
CRs detected by this paper	167844	184019
CRs detected by L.A.COSMIC	163434	173963
efficiency of this paper	73.8%	80.9%
efficiency of L.A.COSMIC	71.9%	76.4%
False detections of this paper	5820	16626
False detections of L.A.COSMIC	559414	38912

Equations

$$P(x)=\alpha e^{-\frac{|x-\beta {|}^{\delta }}{\delta {\gamma }^{\delta }}},$$

$$I_{i,j}^{(2)}=I_{int[(i+1)/2],int[(j+1)/2]},$$

$$\text{L}=\mathcal{L}\ast I^{(2)},$$

$$\mathcal{L}=\frac{1}{4}\left(\begin{array}[]{ccc}0&-1&0\\ -1&4&-1\\ 0&-1&0\end{array}\right).\$$

$${\text{L}}_{i,j}^{+}=\frac{1}{4}({\text{L}}_{2i-1,2j-1}+{\text{L}}_{2i-1,2j}+{\text{L}}_{2i,2j-1}+{\text{L}}_{2i,2j}),$$

$$N_{m5}=\frac{1}{g}\sqrt{g(I_{m5})+{\sigma }_{rd}^2},$$

$$S=\frac{{\text{L}}^{+}}{2N_{m5}}.$$

$$S^{{}^{\prime }}=S-S_{m5}.$$

$$I_{xy}$$

$$\Delta \mathcal{F}$$

$$\begin{array}[]{cccccc}A_{1}&=&\frac{1}{2}\left(\begin{array}[]{ccc}0&0&0\\ -1&0&1\\ 0&0&0\end{array}\right)&A_{2}&=&\frac{1}{2}\left(\begin{array}[]{ccc}0&1&0\\ 0&0&0\\ 0&-1&0\end{array}\right)\\ A_{3}&=&\frac{1}{2\sqrt{2}}\left(\begin{array}[]{ccc}-1&0&0\\ 0&0&0\\ 0&0&1\end{array}\right)&A_{4}&=&\frac{1}{2\sqrt{2}}\left(\begin{array}[]{ccc}0&0&1\\ 0&0&0\\ -1&0&0\end{array}\right)\end{array},$$

$$G=\frac{1}{2}(|A_1\ast \mathcal{F}|+|A_2\ast \mathcal{F}|+|A_3\ast \mathcal{F}|+|A_4\ast \mathcal{F}|).$$

$$N_1=\frac{1}{g}\sqrt{g\mathcal{F}+{\sigma }_{rd}^2},$$

$$N_2=\frac{1}{g}\sqrt{g\mathcal{F}+{\sigma }_{rd}^2}+G.$$

$$D_1=\frac{(I-\mathcal{F})}{N_1},$$

$$D_2=\frac{(I-\mathcal{F})}{N_2}.$$

$$B=\left(\begin{array}[]{ccc}1&1&1\\ 1&0&1\\ 1&1&1\end{array}\right)\ ,$$

$$D_m=[(I-\mathcal{F})(1-M{}^{\prime }){]}_{m3},$$

$$D_1^{{}^{\prime }}=\frac{(I-\mathcal{F}-D_m)}{N_1},$$

$$D_2^{{}^{\prime }}=\frac{(I-\mathcal{F}-D_m)}{N_2}.$$

$$\frac{I-I_{m5}}{N_{m5}}>3,$$

$$\varphi =\frac{f_{CR}}{\sqrt{f_{clean}+{\sigma }_{rd}^2}},$$

Full text

techniques: image processing

1 Introduction

Cosmic rays (CRs, hereafter) are high-energy particles that generate randomly distributed, large signals on charge-coupled devices (CCDs), which could affect the measured fluxes of astronomical objects if not detected or removed properly. Generally, CRs are removed by combining three or more exposures of the same field (Windhorst, 1994; Freudling, 1995; Fruchter, 1997; Gruen, 2014; Desai, 2016) , as they are unlikely to hit the same pixel in more than one exposure. However, multiple exposures are not always available. Furthermore, there are certain situations in which CR detection in single exposures is desired, such as in time domain studies.

Various methods have been developed for identifying and replacing CRs in CCD data of single exposures, including median filtering (e.g., Dickinson’s IRAF tasks QZAP, XZAP, and XNZAP), applying a threshold on the contrast (e.g., IRAF task COSMIC-RAYS ), trainable classification (Murtagh, 1992; Salzberg, 1995; Bertin, 1996), convolution with adapted point-spread functions (PSFs; Rhoads, 2000), Laplacian edge detection (van Dokkum, 2001), analysis of the flux histogram (Pych, 2004) and a fuzzy logic-based method (Shamir, 2005). All of the median filtering or PSF methods remove small CRs from well-sampled data effectively, but problems arise when CRs affect more than half the size of the filter or when the PSF is smaller than the filter (van Dokkum, 2001). All of the methods listed above are designed for photometric data except those of van Dokkum (2001) and Pych (2004), which work for long-slit spectroscopic data.

Farage (2005) made a comparison between different methods including the IRAF script JCRREJ2 of Rhoads (2000), the IRAF routine L.A.COSMIC of van Dokkum (2001), the C script of Pych (2004) and the IRAF task XZAP on photometric images. In that paper, Farage concluded that L.A.COSMIC provided the best performance, with a detection efficiency of 86% on the real data sample, whereas other methods could at most detect 78% of the CRs. Increasing object density reduces the efficiency of detection (Farage, 2005), which is unfortunately unavoidable in multi-fiber spectroscopic data where the signals are always dense. Although L.A.COSMIC efficiently detects CRs, it replaces the identified CR candidates with the median value of the surrounding good pixels (van Dokkum, 2001), which is improper when the CR hits are on the ridge or slope of the profile.

There have been no specific efforts to solve this problem on multi-fiber spectroscopic data, which present distinct challenges compared to photometric data. Multi-fiber images do not have clear isolated point or extended sources as in the photometric data, and the long stripe-like multi-fiber spectra occupy large contiguous regions so that the available area for the local “background” is much smaller than in the photometric data. Methods with median filtering or interpolation of neighboring pixels are less effective in this case.

We present an algorithm to detect and replace CRs for Large Sky Area Multi-Object fiber Spectroscopic Telescope (LAMOST) single-exposure images based on a two-dimen- sional (2D) profile fitting of the spectral aperture. We first pick out CR candidates with Laplacian edge detection and construct a 2D function to fit the image profile in small segments along the spectral trace with these candidates masked out; the final CR list is generated by comparing the fitting residual with a noise model depending on position, and the CR polluted pixels are replaced with the corresponding value of the 2D function. This method is applied to the data processing of LAMOST; in principle, it can also be used for other multi-fiber spectral data after minor modification.

We describe LAMOST data in §2. The algorithm is explained in §3. In §4, we give some examples and analyze the properties of the algorithm. Finally, in §5 we summarize our work.

2 LAMOST Data

LAMOST (Cui, 2012) is a fiber spectroscopic telescope equipped with 4000 fibers feeding 16 spectrographs. Each spectrograph, holding 250 fibers, is split into blue (3700–5900 Å) and red (5700–9000 Å) arms by a dichroic mirror. Groups of 250 spectra are recorded by two $4k{\times}4k$ CCDs at the blue and red end, respectively. The typical duration for a single LAMOST exposure ranges from 600 to 1800s, depending on the target brightness and weather conditions. A considerable number of CRs hit the images during the exposure; for example, in a typical 1800s image, the number of pixels polluted by CRs is about $2\times 10^{4}$.

The size of a LAMOST image, of which the dispersion direction is along the vertical direction, is 4096$\times$4136 pixels. In the spatial direction, the typical distance between two adjacent fibers is 15$\sim$16 pixels. As shown in Figure 1, the cross section of the fiber profile in the spatial direction could be well described by a Sérsic function (Sérsic, 1968; Clewley, 2002):

[TABLE]

where $\alpha$, $\beta$, $\gamma$, $\delta$ are parameters to be derived. The typical full width at half maximum (FWHM) of the profile is about 7$\sim$8 pixels. If $d=|x-\beta|$ is the distance from a given pixel to the fiber profile center in the row (or horizontal/spatial) direction, according to Figure  1, the flux at $d=8$ is less than 0.01% of those at the profile peak. To avoid fiber to fiber cross talk, the magnitude range of objects observed in one LAMOST observation is constrained to be less than 5 magnitudes. In the extreme case, the contribution from the 5 magnitude brighter neighbor to the pixel at $d=8$ could be ignored, so the fluxes in the pixels of $d\leq 8$ could be considered as the flux from the fiber itself. In 2D data reduction, $d=8$ is chosen as the aperture for spectrum extraction. The spectral resolution of LAMOST is about 1800, which corresponds to a FWHM about 5 pixels in the dispersion direction. The PSF changes gradually with position on the CCD chip, but could be considered as constant in a small region (e.g., 20 pixels); we will take advantage of this characteristic to improve the cosmic ray rejection.

3 Cosmic Ray Detection and Rejection

CRs are detected and replaced in three steps. First, we use Laplacian edge detection (van Dokkum, 2001) to generate a raw CR candidate list. Second, for each fiber, pixels within $d=8$ of the fiber trace center are divided into small blocks; each block is then fitted by a 2D profile with those raw CR candidates masked out. The final CR list is determined by comparing the fitting residual with a noise model considering both the intrinsic noise and the uncertainty introduced by profile fitting; the pixels polluted by CRs are replaced by the corresponding fitted value. The details are as follows.

3.1 Laplacian edge detection

Laplacian edge detection has been widely used for highlighting boundaries in processing digital images (e.g., Gonzalez, 1992). van Dokkum (2001) was the first to apply the method to detect CRs in astronomical images. Their publicly available program, L.A.COSMIC, successfully detects CRs in both photometric and long-slit spectroscopic images. We use a similar method to that in Section 3 of van Dokkum (2001) to pick out the raw CR candidates. Since further details can be found in that paper, only basic steps are listed here.

The original image $I$ with the size of ${n_{x}}\times{n_{y}}$ is subsampled into ${2n_{x}}\times{2n_{y}}$:

[TABLE]

where $i=1,\cdots,2n_{x}$ and $j=1,\cdots,2n_{y}$. The subsampled image is then convolved with a Laplacian kernel:

[TABLE]

where $\mathcal{L}$ is the Laplace operator, $*$ denotes convolution. The Laplace operator $\mathcal{L}$ in the above convolution is

[TABLE]

Since CRs are positive in $\L$, all the negative values are set to zero. $\L$ is then resampled to ${n_{x}}\times{n_{y}}$:

[TABLE]

where $i=1,\cdots,n_{x}$ and $j=1,\cdots,n_{y}$.

The original image is median filtered with a ${5}\times{5}$ box to construct the noise model

[TABLE]

where $g$ is the gain in electrons per ADU, $I_{m5}$ is the image median filtered by a ${5}\times{5}$ box, and $\sigma_{rd}$ is the readout noise in electrons. The Laplacian image is then divided by the noise model and the subsampling factor to obtain the deviations from the expected Poisson fluctuations:

[TABLE]

All structures that are smooth on scales of $\geq{5}$ pixels are removed by a ${5}\times{5}$ median filter:

[TABLE]

All pixels that meet $S^{{}^{\prime}}>\sigma_{lim}$ are identified as CR candidates, where $\sigma_{lim}$ is a given threshold. We adopt $\sigma_{lim}=4.5$, similar to L.A.COSMIC.

3.2 Fiber profile fitting

Fiber traces are closely aligned on multi-fiber spectral images. In contrast to photometric images, useful signals are quite fully rather than sparsely distributed. As pointed out by Farage Farage (2005), the increasing object density will certainly reduce the efficiency of the methods designed for photometric data. The ramp on either side of the ridge of the fiber trace is quite steep, so it is hard for the edge detection method to discriminate between real CR hits and good pixels on the ramp, leading to a drop in the detection rate and a rise in the false detection rate. Furthermore, replacing the CR polluted pixels with the median of the surrounding pixels is seemingly unsafe. Along the fiber trace, the shape of the PSF changes slowly. If the PSF is well sampled, then CR discrimination could be improved by its shape difference from the PSF.

For an image $I$ of ${n_{x}}\times{n_{y}}$ pixels with $n_{f}$ fibers, the pixels close to the fiber trace center contribute the most to the extracted flux. CRs hitting on these areas will introduce large errors in the final spectra, while those in the trough between fibers have much less impact. Consider a small spectral segment centered on column $[c_{kj}]$ and row $j$, where $k=1,\cdots,n_{f}$, $j=1,\cdots,n_{y}$ and $[c_{kj}]$ is the trace center of the $k$th fiber at row $j$($[]$ denotes the round off of the quotation). Since both the trace center and the shape of the PSF change slowly inside the segment, the shape of the segment could be fitted with a product of two orthogonal vectors:

[TABLE]

where $S(x)$ is the fiber profile in the spatial direction, $P(y)$ is a polynomical to describe the flux variation in the dispersion direction, and $\varepsilon_{xy}$ is noise. The size of $x$ is the same as the aperture for flux extraction, which is set to $d=8$ for LAMOST, as discussed in Section 2. The size of $y$ is chosen to be small enough to keep $P(y)$ smooth but larger than the size of single CR hits, so that $P(y)$ can be fitted with a low-order polynomial and the CR polluted pixels could be better estimated by interpolation. For LAMOST, the segment size is set $17\times 9$, i.e., $x=[c_{kj}]-8,\cdots,[c_{kj}]+8$ and $y=j-4,\cdots,j+4$. We do not try to fit the pixels in the bottom of the valley between fibers, since they contribute little to the final spectrum.

The shape of $S(x)$ is determined by the output pupil of the fiber and instrument distortion. Although the typical shape of $S(x)$ could be described by a Sérsic function (Eq. 1), the actual shape deviates occasionally from the analytic function when the optical distortion is large at the edge of the image or the coupling between the fiber output pupil and the slit is imperfect. Due to the above reason, $S(x)$ is constructed with an empirical profile rather than an analytic function. All the profiles at $y=j-10\sim{j+10}$ are first normalized, center justified in sub-pixel scale and then averaged to derive $S(x)$ with the CR candidates masked out. Fixing the form of $S(x)$, the polynomial coefficients of $P(y)$ are derived by least-square surface fitting to the flux in the segment with the CR candidates masked out. A fitted image $\mathcal{F}$ is generated after all segments are handled.

3.3 Cosmic ray selection

A new list of CRs is generated by comparing the noise model with the residual of the image fitting without reference to the old CR list. The noise or uncertainty of our method comes from two parts: one is the intrinsic noise of the input signal, i.e., Poisson noise from the object and the readout noise from the CCD circuit; the other part of the noise comes from the defect of the profile fitting, which is larger when the profile changes more dramatically. Basically, the first part is related to time and the second part is related to position, which could be illustrated as:

[TABLE]

The first term is sufficient to pick out CRs for most cases when the fit is good, yet it is necessary to add the position- dependent term to avoid false detections in regions where the fit is not perfect. The following steps are implemented to reject the CRs:

The position dependent terms in 10 are calculated by the average gradients at each pixel. Convolving the fitted image  with the following four arrays:

[TABLE]

a gradient array can be derived as

[TABLE] 2. 2.

Noise models are constructed with and without the second term in Equation 10, respectively,

[TABLE]

and

[TABLE]

The noise-weighted differences between the input image $I$ and the fitted image $\mathcal{F}$ are defined accordingly:

[TABLE]

and

[TABLE]

All pixels with $D_{1}>20$ or $D_{2}>3$ are marked as CR candidates in this step. A mask array $M$ is generated with the CR polluted pixels set to 1. 3. 3.

By setting the previous limits, the number of fake CRs is greatly reduced, but real CRs with a low signal-to-noise ratio (S/N), most of which are indiscernible from noise, are blocked as well. Considering the consecutive pixels occupied by a certain CR hit, the pixels at the edge of the CR hit are more likely to be rejected due to lower signal, though they should have higher probability to be real than the single-pixel event. So a lower limit for those neighbouring pixels will raise the detection rate. To do this, all of the neighboring pixels are first added back to the CR list by convolving the mask array $M$ with

[TABLE]

for the pixels in the expanded CR list, if the corresponding $D_{1}>2$ or $D_{2}>2$, then they will be added to a new mask array $M{{}^{\prime}}$. 4. 4.

Assuming there are sufficient CRs masked out, the residual of the fit $(I-\mathcal{F})$ will represent the actual difference between the original CR-free image and the fitted image. The difference can be added back to compensate for imperfect fitting, and a more accurate fit will help to raise the CR detection rate as follows.

The difference is derived by median filtering the residual $(I-\mathcal{F})$ with the CR candidates $M{{}^{\prime}}$ masked out:

[TABLE]

where m3 denotes a $3\times 3$ median filter. And the new noise-weighted residual arrays will be

[TABLE]

and

[TABLE]

All pixels with $D^{{}^{\prime}}_{1}>10$ or $D^{{}^{\prime}}_{2}>3$ or $D_{1}>20$ or $D_{2}>3$ are masked as CR candidates in this step. 5. 5.

Rerun to Step 3 and confirm the final CR candidates.

We do not try to fit the pixels that are either in bad fibers or in the $d>8$ gaps between fibers because they contribute little to the final extracted spectra. For those pixels, $\mathcal{F}$ is set to 0 and the CR candidates are selected by simply requiring the noise- weighted difference between the original image $I$ and the $5\times{5}$ median filtered image $I_{m5}$ to be larger than 3:

[TABLE]

where $N_{m5}$ is the same as in Equation 6.

Combining the above CR candidates, the final CR mask is generated and the value of each CR polluted pixel is replaced by the corresponding value in the fitted image $\mathcal{F}$.

3.4 Additional features

With the profile fitting method, other bad pixels such as the inherent damaged pixels could be replaced with a reasonable value once an initial bad pixel map is known.

On a Dell Precision T5500 (eight 2.0 GHz CPUs), the IDL implementation with one single-threaded processor requires about 20 minutes for an image of ${4096}\times{4136}$ pixels. Most of the time is spent on the 2D profile fitting and the timescales linearly with the number of fibers and image size in the dispersion direction. As the current version of our program is not parallelized, the execution time on the current computer is equivalent to that of a single-core processor and could be greatly reduced after software parallelization.

4 Examples and Application

In this section, tests with both simulated data and real data are carried out to illustrate the performance of our method. In these tests, our primary concerns are the following factors: how many pixels of real CRs are detected (efficiency), how many pixels are falsely detected as CRs (false detection rate), and the accuracy of the CR replacement. A better method should have higher efficiency, lower false detection rate, and proper restoration of the pixels polluted by CRs.

4.1 Artificial Images

For multi-fiber spectral observations, to increase the observation efficiency as well as to avoid fiber to fiber cross talk, it is usually a good strategy to divide the targets into different plates according to their brightness so that the S/N of the targets in the same plate are similar in the same exposure time. For plates with bright magnitude, the exposure time is short, so the strength of the sky spectrum is low, but the S/N of the objects is high. On the contrary, the exposure times for the faint plates are long, so the sky spectrum is strong and the object spectrum is weak. Since the efficiency strongly depends on the brightness contrast between the CR and the target, plates with different target brightness are simulated to test our method under different situations.

First, we generate a pure CR image of ${4136}\times{4096}$ pixels with 20,000 CR hits (approximately 10 times of those in a 30 minutes exposure LAMOST image). The shape of each CR is set to be an ellipse with the major axis randomly distributed in 1$\sim$10 pixels and the minor axis ranges from 1 to 3 pixels, by which almost all kinds of CRs in the real image can be simulated. The direction of the major axis is randomly distributed in 0$\sim 360^{\circ}$ and the intensities are uniformly distributed between 0$\sim$20000 ADUs. All CRs with flux less than 5, which is at the readout noise level, are set to 0. In total, 227,451 pixels are polluted by CR hits.

Second, to study the method’s performance with different target brightness, two CR-free images are generated by combining three consecutive LAMOST exposures of the same targets. The first image (IMG600) is combined from three 600 second exposures in which the spectra are dominated by the strong smooth continuum from the bright objects and the sky emission lines are relatively weak. In the second one (IMG1800), the sky emission lines are more prominent (due to a longer exposure time of 1800 seconds) and the object continuum is relatively weak (due the faint magnitude). The final test images are generated by adding the pure CR image to the CR-free images.

Both IMG600 and IMG1800 are tested by our method and the IDL version of L.A.COSMIC, respectively. The results are summarized in Table 1. For IMG600, the efficiency of our method is 1.9% higher (73.8% versus 71.9%) and the number of false detections is 5820, two orders of magnitude lower than that of L.A.COSMIC. For IMG1800, our efficiency is 4.5higher(80.9% vies 76.4%) and the number of false detections (16,626) is less than half that of L.A.COSMIC(38,912).

The efficiency of both methods rises more than 4% from IMG600 to IMG1800. The reason is that the efficiency, for those CRs falling coincidentally into the same pixel with the object spectrum, decreases with the increasing photo noise, while in this case, the CRs are the same in both simulations but the object spectra are much brighter in IMG600 therefore the noise is larger in IMG600 than in IMG1800. Compared with our method, L.A.COSMIC is prone to mistake the wings of the bright profiles as the sharp edges of CRs, especially when the contrast between the background and the profile peak is high, as in IMG600. In IMG1800, as the object brightness decreases, the contrast and thus the number of false detections drops. Our profile fitting method successfully bypasses this sharp edge trap in IMG600, reducing the huge number of fake detections to a reasonable level. As shown in Table 1, for our method, the number of false detections doubles from IMG600 to IMG1800. There are two reasons for this problem. First, as the SNR of the spectrum becomes lower, more faint pixels are mistaken as CRs (as can be seen from Figure 4). Second, as the exposure time increases, the intensity of the sky emission lines increases, but the intensity of the underlying object spectra decreases (for the magnitude gets much fainter). In this case, the relative change at some of the exponential wings of the strong sky emission lines becomes too dramatic to have a good polynomial fit; the larger residual induced by the improper fitting leads to an increase of false detections.

If we denote the SNR of a CR polluted pixel as

[TABLE]

where $f_{CR}$ and $f_{clean}$ are the fluxes from the pure CR and the CR-free image, respectively, $\sigma_{rd}$ is the readout noise, then Figure 2 shows the detection efficiency against $\phi$. Most of the undetected CRs are those with low $\phi$. The efficiency remains high for $\phi>10$ then drop quickly when $\phi<10$. The efficiency of our method is higher than that of L.A.COSMIC  in all situations except for $\phi<2$, where the CRs are too weak to be separated.

In Figure 3, the recovered fluxes of IMG600 are compared to the corresponding fluxes of the CR-free image to see how well our CR replacement works; as shown in the left and the right panel, almost all the replacements properly follow the original fluxes. Also shown in the middle panel is the replacement performance of L.A.COSMIC; most of the replacements are good, but the scatter is larger especially in the high flux region, which is not unexpected, since its replacement method is not specially designed for multi-fiber spectra. The performance of both methods on IMG1800 is similar to IMG600.

Since the replacement of the false detections also changes the flux, causing errors in the spectrum, it’s necessary to test the replacement on those falsely detected CRs. L.A.COSMIC produces too many falsely detected CRs to be comparable with our method in IMG600, so only the results of IMG1800 for both methods are shown in Figure 4. As shown in the picture, the results of both methods deviate from the true value. Though our method systematically underestimates the flux, most of our replaced fluxes concentrate within 80% of their true value and the true fluxes of most of the pixels are low, so the influence on the extracted spectrum should be small. L.A.COSMIC results show a large variation, with a large number of pixels shifting from the true value to very low fluxes. Some examples are demonstrated in Figure 5.

For spectroscopic data, the extracted spectra are more important than the flux of individual pixels on the 2D image. Figure 6 compares the fluxes of the CR polluted part of the extracted spectra with the CR-free spectra in different situations: from top to bottom, it shows the falsely detected, the properly detected, and the undetected CRs, respectively; in all cases, the average difference between the CR corrected spectra and the original spectra is less than 2.2%, as shown from the distribution in the right column of Figure 6 . Figure 7 shows an extracted spectrum sample; the residual of the CR correction is within a few percent.

4.2 Real Data

We test our algorithm with real data from the LAMOST survey. Figure 8 shows a part of an 1800s LAMOST image (left panel) and its reconstruction by our method (right panel). Visual inspection of our reconstructed image shows that most of the CR hits are properly removed. For a further comparison, the extracted flux of the CR detected pixels vs the flux of the classical multi-exposure-combination method are demonstrated in Figure  9. The results are comparable to the simulations in Figure 6, except that the scatter is a bit larger. The reasons for the larger scatter could be the following: first, the simulated image in Figure 6 has a higher SNR than the real data; second, the falsely detected CRs cannot be discriminated from the true CRs in this test, so the scatter should be larger than the true-CR-only situation; and third, the sky flux varies between exposures and the object flux gathered by LAMOST varies with telescope pointing, which make the combined image deviate from individual exposures, leading to a larger scatter.

5 Conclusion

We present a method for detecting and removing 2D profile fitting to each segment. A new cosmic ray list is generated by comparing the fitting residual with a noise model depending on both the intrinsic shot noise and the relative position in the profile. We finally produce a more accurate cosmic ray mask table and more reasonable substitution values for CR polluted pixels. The method is tested by both simulations and real data; the results show that our method has a high detection rate, low false detection rate, and proper replacement of the CR polluted pixels.

Since this method fits the 2D profiles of the fiber spectro- scopic data, which are different from the photometric PSF, it cannot be applied to photometric data. However, it can be used in slit spectroscopic data after minor modifications. The code and samples are available at http://lamostss.bao.ac.cn/~bai/crr.

Z. Bai acknowledges the support of the National Natural Science Foundation of China (NSFC) (grant no. 11503054). H. Zhang acknowledges the support of NSFC Key Program (grant no. 11333004) and the National Key Basic Research Program of China (grant 2014CB845700). The Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectro- scopic Telescope, LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.