Optimal Cosmic-Ray Detection for Nondestructive Read Ramps
aa r X i v : . [ a s t r o - ph . I M ] S e p Draft version November 13, 2018
Preprint typeset using L A TEX style emulateapj v. 8/13/10
OPTIMAL COSMIC-RAY DETECTION FOR NONDESTRUCTIVE READ RAMPS
Rachel E. Anderson and Karl D. Gordon
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218 (Dated: November 13, 2018)
Draft version November 13, 2018
ABSTRACTCosmic rays are a known problem in astronomy, causing both loss of data and data inaccuracy.The problem becomes even more extreme when considering data from a high-radiation environment,such as in orbit around Earth or outside the Earth’s magnetic field altogether, unprotected, as willbe the case for the
James Webb Space Telescope ( JWST ). For
JWST , all the instruments employnondestructive readout schemes. The most common of these will be “up the ramp” sampling, wherethe detector is read out regularly during the ramp. We study three methods to correct for cosmic raysin these ramps: a two-point difference method, a deviation from the fit method, and a y-interceptmethod. We apply these methods to simulated nondestructive read ramps with single-sample groups,and varying combinations of flux, number of samples, number of cosmic rays, cosmic-ray location inthe exposure, and cosmic-ray strength. We show that the y-intercept method is the optimal detectionmethod in the read-noise-dominated regime, while both the y-intercept method and the two-pointdifference method are best in the photon-noise-dominated regime, with the latter requiring fewercomputations.
Subject headings:
Data analysis and Techniques, Astrophysical Data, Astronomical Techniques INTRODUCTIONAdvances in technology that allow us to observe fainterobjects, build more complex systems, and send telescopesfurther into space have challenged us to continue to im-prove our calibrations. This includes detection methodsfor cosmic rays (CRs) or any charged particle that addsa jump to the data.
JWST ’s orbit at the second Earth-Sun Lagrange point (L2), which allows passive coolingof the telescope to ∼
50 K, puts it outside the protec-tive mantel of the Earth’s magnetic field. This couldmake CRs a serious problem. Furthermore, newer in-frared telescopes will have a lower read-noise, and thuswe will detect lower CRs. Finally, long observing timeswill be necessary to complete many of the scientific goalsof the
JWST , causing CRs to be an even larger problem.From a study by Robberto (2009a) we can calculate thatin every 2000 s, on average, 13% of the pixels on the
JWST
HgCdTe detectors and 25% of the pixels on the
JWST
Si:As detectors will be affected by CRs. Thesevalues could be even greater, since this study does nottake into account secondary particles. For comparison,onboard measurements of the
NICMOS camera on the
HST show that about 10% of the pixels show a CR hitfor every 2000 s of integration (Viana and Wiklind et al.2009). These CRs will include low-energy CRs from sec-ondary particles and the Sun, and higher-energy galacticCRs (Robberto 2009c). Clearly, a reliable method to de-tect both low- and high-energy CRs is needed.In this article we will discuss CR detection methodsfor infrared data, which use nondestructive read ramps.Offenberg et al. (2001) found that in order to correct forCRs, the nondestructive readout scheme was most effi-cient. By integrating the charge on each pixel in thisway it is possible to calculate the slope of these ramps incounts per time to get the flux of the sky (Rieke 2007). [email protected], [email protected]
We refer to the integrated charge as a ramp made upof a specified number of samples. We discuss the slopeof these ramps as the calculation of the flux. Finally, aCR affects the ramp between samples, however we use‘sample number’ to refer to the first sample after the CRhit.Correcting for CRs in ramps is not a new prob-lem (Offenberg et al. 1999). However, the advent of largenondestructive arrays in orbit coupled with ground-basedprocessing of the ramps means, when it comes to calcu-lating the slope of ramps, that there are more optionsavailable for CR detection. In addition to CRs, noise isadded to the data by the detector and readout electron-ics as well (Fanson 1998, Tian et al. 1996, Rieke 2007).In this article we focus on photon-noise and read-noisefor our correlated and uncorrelated noise, respectively;however, we are not restricted to just these two terms.The work in this article is general and can account for anycorrelated and uncorrelated noise sources. Therefore, ourquestion is, What is the best we can do at finding CRsin a ramp, given the noise in the ramp?To understand and test various CR detection methodsfor ramps, we simulate infrared data as described in Sec-tion 2. Anytime a slope is calculated for the precedingprocess or for the CR detection methods, we do so usinglinear regression for data with correlated and randomuncertainties (in consideration of the photon-noise andthe read-noise, respectively), as described in Section 3.In Section 4 we propose three techniques to detect CRs.The first is a two-point difference method, the second isa deviation from the fit method, and the third is a y-intercept comparison. There are various conditions thatcan hinder/aid in CR rejection (e.g., number of samples,slope, number of CRs, size of CRs, and location of CRs);therefore, we aim to study combinations of these and tofind which algorithm behaves the best under differentconditions. These results are presented in Section 5. A Anderson & Gordon
Parameters ValuesSlope 70.0 e − /sY-Intercept 21000.0 e − Number of Samples 40Sample Time 27.7 sRead Noise 16.0/ √ e − /sample Table 1
Parameters Used to Create MIRI Ramps. discussion on our findings is in Section 6. Concludingremarks can be found in Section 7. SIMULATING NONDESTRUCTIVE READRAMPSIn order to test how well a CR detection method works,we have to test the method on a known CR. For thatreason, we have simulated ramps with photon-noise andread-noise added in, to which we can either add CRswith known location in the exposure and known am-plitudes or leave the ramps clean (CR-free). We buildthese ramps with the slope, y-intercept, number of sam-ples, time between samples (sample time), and read-noiseas inputs, using parameters given in Table 1 as guide-lines (slope and number of samples will be changed inthis study). The only value that is instrument-specificis the read-noise, which was chosen to match the ex-pected value of the Mid-Infrared Instrument (
MIRI ) onthe
JWST . Note that we do not group any of the sam-ples by taking a weighted average or coadding. Further-more, we assume uniform sampling (i.e., constant timebetween samples), as that is what is used by the
JWST and also so that calculations are more efficient. We alsoassume that a nonlinearity correction has already beenperformed with no error, and we do not account for theeffects of quantization noise. Although we are simulat-ing
MIRI detector parameters, this is only an example.The methods discussed in this article will apply to anyother nondestructive read data, including those for the
JWST near-infrared detectors (
Near-Infrared Camera [ NIRCam ] and
Near-Infrared Spectrograph [ NIRSpec ])the
HST infrared detectors (
Wide Field Camera 3 - IR [ WFC3-IR ] and
Near-Infrared Camera and Multi-ObjectSpectrometer [ NICMOS ]) and ground-based detectors,by simply changing the values in Table 1. If the datayou would like to simulate do included grouped samples,however, some revisions will have to be made.At time t = 0 the counts y are equal to the y-intercept. This is not considered part of the ramp, be-cause t = 0 is when the reset occurred, but we startreading at time equal to one sample time, t .Each signal is comprised of the counts from the pre-vious signal plus the additional signal, the photon-noisefrom the additional signal, and the read-noise. To con-struct the ramps, we follow this procedure:1. For each ramp calculate the expected signal (i.e.,no noise, clean), s e , defined as s e ≡ mt s , (1) Although we use the value for the y-intercept given in Table 1,this is an arbitrary value and does not change our calculations. where m is the input slope, and t s is the sampletime.2. For each sample calculate the sum of the expectedsignal and the unique photon-noise of this signal( s e + σ n p i ) from a Poisson distribution with λ = s e .Add this to the signal from the previous sample (ory-intercept if it is the first sample). This step isshown in equation 2. Since we have only added thecorrelated noise, and we still need to add the uncor-related noise, we have called the signal y ∗ i insteadof y i ). y ∗ i = y ∗ i − + s e + σ n p i . (2)3. If there are to be one or more CRs to the ramp,simply add the electrons with the expected signal: y ∗ i = y ∗ i − + s e + σ n p i + CR mag , (3)where CR mag is the magnitude of the CR.4. Read noise, denoted by n r , is due to readout elec-tronics; consequently, it is uncorrelated. Therefore,when all samples are populated, add a unique read-noise, σ n r i , to each sample (equation 4): y i = y ∗ i + σ n r i . (4)For each sample σ n r i is taken from a Gaussian with n r as the standard deviation.The uncertainties for the samples in the ramp are cal-culated by adding the photon-noise and read-noise inquadrature: σ y = q n r + n p , (5)where n p is the Poisson noise of the expected signal in theramp; thus, n p = √ mt s , and n r is taken from Table 1.We have chosen to give each sample equal weighting,rather than choosing the photon-noise to be the Poissonnoise of the counts in the ramp, so that we can improveour results by avoiding a bias based on sample number. LINEAR FIT ALGORITHM FOR DATA WITHCORRELATED ERRORSFor two of the three CR detection methods, it isvery important that we use the best calculation of theslope and y-intercept. Therefore we must use all in-formation including correlated and uncorrelated errors.Gordon et al. (2005) explains how to take into accountboth correlated and uncorrelated errors for the slope andy-intercept uncertainties, and here we describe how toinclude this information for the fit as well, using matrixnotation and specifically making use of the covariancematrix. For a refresher on calculating a linear fit usingmatrices as well as the use of the covariance matrix, seeHogg et al. (2010). Every time a slope is calculated inthe CR detection methods described in this article, wedo so using this method.The real power of using matrix notation comes fromthe covariance matrix, C , as it includes informationabout the degree of correlation between samples. Fur-thermore, the slope and y-intercept uncertainties result-ing from using this covariance matrix automatically takeetecting Cosmic Rays in Infrared Data 3into account the correlated and uncorrelated errors, pro-vided they are included in the covariance matrix. C isdefined as: C = σ y c , ... c ,N c , σ y ... c ,N ... ... ... ...c N, c N, ... σ y N . (6)The values on the diagonal of C are the uncertain-ties of the y i , whereas the off-diagonal elements are thecovariance between the different y i . If the data are un-correlated, then the off-diagonal matrix elements are allzero, c i,j = 0 , i = j . However, we do have correlateddata due to the photon-noise.When constructing C for correlated errors, we followthe logic of Fixsen et al. (2000) and think of C as thesum of two matrices: one for the photon-noise, P , andone for the read-noise, R . Fixsen et al. (2000) definethese matrices for the case where the data are the two-point difference of the samples (the photon-noise is notcorrelated, and the read-noise is correlated). However,for fitting a line to the individual samples, it is just theopposite. The read-noise is uncorrelated, and thus R isjust n r on the diagonal and zero elsewhere.The diagonal of P is n p , just as it would be for uncor-related data, but since the photon-noise is correlated, theoff-diagonal elements are not zero. An estimate of all ofthe photon-noise in y i is the photon-noise added to eachsample multiplied by i (i.e., in p ). The correlation be-tween samples y i and y j is the photon-noise that is in y i that is also in y j . The elements in P that represent thiscorrelation are p i,j and p j,i . Consider y and y , wherethe photon-noise they share would be the photon-noisein y (2 n p ), plus the photon-noise in the first read, p , .Thus, p i,j = p j,i = kn p + p , , (7)where k = n j, j < ii, i < j . (8)If p , is included in the background, then it is set tozero. We can estimate n p using the slope as we did inSection 2: n p = √ s e . This initial calculation of theslope for s e (eq. 1) is done before any CR correction andwithout taking into account correlated errors.This technique accounts for the diagonal elements aswell. Therefore, c i,j = n p k + p , + n r i = j,c i,j = n p k + p , i = j, (9)where n r is the read-noise.An equation for the noise in a ramp has previously beenderived in nonmatrix form in Rauscher et al. 2007; foran independent derivation and the correct final formula,see (Robberto 2009b). Rauscher et al. formula (eq. (1)in that article, see erratum) has the benefit that it takesinto account grouped samples. However, the benefit ofusing matrix notation with a covariance matrix is thatwe can add other noise terms easily by calculating theappropriate covariance matrix (like we did for P and R ) and sum all to get C . It would be more difficult to addother noise terms to the Rauscher et al. (2007) equation,and like our method here, they only include read-noiseand photon-noise.3.1. Validating with Simulations
To demonstrate how well this calculation of the un-certainties fits simulated data, we followed the idea fromFigure 16 in Gordon et al. (2005) and simulated 10,000ramps each with read-noise only, photon-noise only, andboth, then we used the covariance matrix to calculatethe y-intercept and slope uncertainties. This is given inFigure 1. The dashed and dotted lines are the uncertain-ties with read-noise and photon-noise only, respectively.The solid lines are the uncertainties in the ramps whereboth read-noise and photon-noise were added. Noticethat you can see where the transition is between the read-noise-dominated regime and the photon-noise-dominatedregime using these plots. The circles are the uncertain-ties calculated using the method described in Section 3,which fit the data perfectly.In addition to using a covariance matrix, other popularmethods to fitting a line to ramps have used uncertain-ties as the weights (1/ σ ), equal weights, or optimumweights, as discussed in Fixsen et al. (2000). To demon-strate how these methods compare, we simulated 10,000ramps, each with the same input slope (before the noisewas added) and then tried to retrieve these slopes usingthe various weighting schemes (Figure 2). The standarddeviation in the slope fits is plotted in Figure 2 (a). Wesee that the covariance matrix and optimal weighting re-sult in a higher signal-to-noise ratio (S/N) than the othertwo methods. Optimal weighting estimates a covariancematrix, so it is expected that they show similar results.In Figure 2 (b) we show the ratio of the average cal-culated slope to the input slope, with unity subtracted.Note that the curve in this plot is nonsmooth due tothe finite number of trials. The error on the slope cal-culation is very small, with the calculated slope at most0.004% from the input slope. For various slopes we madea histogram of the distribution of these ratios and foundthat they were symmetric and centered on the average.Therefore, optimal weighting and using a covariance ma-trix with correlated and uncorrelated errors produces alower standard deviation of the results, while all methodsproduce a similar slope estimate with 10,000 trials. CR DETECTION METHODSAll three CR detection methods have a unique algo-rithm. However, all follow the same procedure to cal-culate a slope once the CRs are found, which we nowoutline:1. Detect CRs, always looking for the largest outliersfirst.2. Calculate the slope for the resulting ramp seg-ments before and after the CR events (we will referto these as “semiramps” from here on, followingRobberto 2008).3. Calculate the final slope of the entire ramp. If thereis one CR or more, do this by taking the weightedaverage of the slopes of the semiramps (see also thediscussion in Section 6). Anderson & Gordon -4 -3 -2 -1 Slope [e-/s]10 -1 Y - I n t e r c e p t U n c . [ e - ] n r to n p transition (a)Read Noise, n r Photon Noise, n p BothBoth (Matrix C) 10 -4 -3 -2 -1 Slope [e-/s]10 -3 -2 -1 S l o p e U n c . [ e - / s ] n r to n p transition (b)Read Noise, n r Photon Noise, n p BothBoth (Matrix C)
Figure 1.
Calculation of the uncertainties for the y-intercept (a) and the slope (b) using the method described in Section 3. Following theidea from Figure 16 in Gordon et al. (2005), for 10,000 each, we generated ramps with read-noise only, photon-noise only, and then both.The uncertainties in the y-intercept and slope given by the data are plotted by the dashed, dotted, and solid lines. Then using correlatedand uncorrelated components in the covariance matrix (e.g., including both read-noise and photon-noise) we calculated the correspondinguncertainties shown by the circles. These calculations match the data perfectly. Furthermore, this plot also shows where the transition isbetween the read-noise-dominated regime and the photon-noise-dominated regime. S t d e v i n S l o p e F i t s (a) 1/ (cid:0) WeightsEqual WeightsOptimum WeightsCovariance Matrix0 200 400 600 800 1000 1200 1400Slope [e-/s] (cid:1) (cid:2) (cid:3) (cid:4) (cid:5) S l o p e c a l c / S l o p e i n p u t - (b) Figure 2. (a) Standard deviation of slope fits versus input slope,for various weighting schemes to point out that using the covari-ance matrix has lower noise (though it is very similar to optimumweighting, as optimum weighting is an estimate of the covariancematrix). (b) Calculated slope to input slope ratio minus one forthe same methods. All methods show the same bias; therefore, thedifference lies only in the S/N shown in (a).
Furthermore, for all methods we use the absolute dif-ference so that we remove outliers in both directions anddo not bias the data (Offenberg et al. 1999, 2001). How-ever, if larger rejection thresholds are used, and thereforethere are no longer as many outliers (i.e., only picking upthe strongest CRs), then one-sided clipping would workjust as well (Windhorst et al. 1994).4.1.
Two-Point Difference Method
With the nondestructive readout method, algorithmsfor CR detection include computing the two-point differ-ences. Offenberg et al. (1999) and Fixsen et al. (2000)describe versions of this method, though both use onlyone rejection threshold for all cases (5 σ and 4.5 σ , re-spectively). Another version was used in the MultibandImaging Photometer ( MIPS ) data reduction algorithmfor
Spitzer (Gordon et al. 2005). We describe a variantof this technique here. C o un t s [ e - ] (a)Data2-Point Diff, d i | d i (cid:6)(cid:7) d | / (cid:8) d (b) Figure 3. method. (a) Solid line and points are thedata, y i , and the dotted line is the two-point differences, d i , of thisdata. (b) Dashed line is the ratio in equation 10 which will becompared with the rejection threshold, r t . The vertical line showsthe interval where the CR hit. For the two-point difference method (hereafter ), we calculate the two-point differences between thecounts in each set of adjacent samples. The largest out-lier is flagged as a CR, given that it fulfills the rejectioncriteria: | d i − µ d | σ d > r t , (10)where d i is the difference between the science data y i +1 and y i , µ d is the median of d , σ d is the uncertainty of d , and r t is the rejection threshold. The median is usedbecause it is more robust than the mean when there areoutliers in the data (Press et al. 1986; Offenberg et al.1999). The median can become a problem if there isquantization noise (i.e., if the slope of the ramp were 0.07 e − /s), however we do not take into account quantizationnoise in this article. Once a CR is identified, the d i thatincludes the CR is removed, and the process is repeatedon the remaining d i until no more CRs are found. Thismethod is depicted in Figure 3.etecting Cosmic Rays in Infrared Data 5When calculating the uncertainty in d , σ d , the mostobvious solution would be to use the standard deviation.However we have found that this does not work well forramps with a low number of samples ( ∼ σ d by using the photon-noise and read-noise addedin quadrature. The photon-noise can be calculated asPoisson noise, but since we are dealing with the two-point difference we use the charge accumulated since thelast sample, rather than the total charge in a sample, sowe use √ s e . Since a CR can contaminate a slope calcula-tion, we can use µ d to estimate s e . Therefore n p = √ µ d ,and σ d = q n p + 2 n r , (11)where the factor 2 is due to the read-noise from each ofthe two samples.4.2. Deviation from the Fit Method
The method we will refer to as the deviation from thefit method (hereafter dev from fit ), is the one usedby
NICMOS (Dahlen et al. 2008) and
WFC3 (Dressel2010). To use this method we fit a line to the ramp usinga covariance matrix as described in Section 3. Then, foreach sample, we calculate the difference to the fit as aratio to the uncertainty in the counts: dev i = y i − f i σ y , (12)where f i is the fit at sample y i , and σ y is the uncertaintyin each sample, defined in equation 5.We then take the first difference of these ratios, andlook for the largest. If it satisfies the criteria: | dev i +1 − dev i | > r t , (13)then y i +1 is flagged as a CR. The ramp is then splitinto semiramps, and this method is applied again to theresulting semiramps.The dev from fit method is illustrated in Figure 4.Note that the background level is not at zero as it wasfor the method, but ∼
10. The reason forthe nonzero background is that while the median (whichwould exclude the CR) is subtracted in the method, in the dev from fit method the fit is sub-tracted, which includes the CR in its calculation. Fur-thermore, in the presence of a CR, the slope will be over-estimated and, therefore so will the photon-noise. Fi-nally, the peak is at a ratio of ∼
60, whereas both the and y-int methods (as we will see) peak atratios of ∼
80. 4.3.
Y-Intercept Method
The idea for the y-intercept method (hereafter y-int ) comes from the
MIPS data reduction algo-rithm (Gordon et al. 2005). The details we describe hereare a variant of that method. For the y-int method,we step through each sample and assume that there isa CR there (the first sample is skipped). We fit a lineto the semiramps before and after the sample with theassumed CR, using a covariance matrix as described inSection 3. At each sample, the x-axis is shifted so thatthe y-intercept is located at the sample number of the as-sumed CR. The only exception is when we assume that C o un t s [ e - ] (a)DataLinear Fit, f i (cid:9) (cid:10) D e v . f r o m F i t (b)dev i |dev i (cid:11) dev i (cid:12) | Figure 4. dev from fit method. The solid line and points arethe data, y i and the dotted line is the linear fit made up of points, f i . In the lower panel the dotted line is the ratio of the differencebetween the fit and the data to the uncertainty, dev i . The dashedline is the two-point differences of the dev i . This is what will becompared with the rejection threshold, r t . The vertical line showsthe interval where the CR hit. there is a CR in the second sample. Since we cannotcalculate a y-intercept for only one point, we shift thex-axis to the first sample, and use the counts in thatsample as the y-intercept, and take the read-noise as theuncertainty. We do the mirror to this when assumingthat there is a CR in the last sample. Then, we takethe ratio of the absolute differences between the two y-intercepts ( b and b ) to the expected uncertainty, σ b (equation 14). | b − b | σ b > r t (14)After stepping through every sample, we look at the sam-ple with the largest ratio. If this ratio is larger than agiven r t , then we flag that sample as a CR. The processis repeated on the resulting semiramps. The method isdepicted in Figure 5.The expected uncertainty, σ b , is calculated from they-intercept uncertainty from the read-noise and the y-intercept uncertainty from the photon-noise added inquadrature. There are two read-noise components, n r and n r , one from the ramp before the assumed CR, andone from the ramp after. These read-noise terms arecalculated by setting the covariance matrix, C , to theread-noise matrix, R , and recalculating the y-interceptuncertainty. The photon-noise component is calculatedas n p = √ s e . Since in the y-int method we are deal-ing with two ramps, we take the weighted average of theslopes to get s e , just as we would if there were actuallya CR found at the assumed location. Therefore, σ b = q n p + n r + n r . (15)In Figure 5 (c) the dashed line representing equation 14is hovering around zero before the CR detection, andthen afterward it increases to ∼
10. This effect comesfrom the difference between the y-intercepts. It is causedby the fact that we shift the x-axis to the sample with theassumed CR, whereas if we shift it to the sample beforethe CR (as we do for a CR in the second sample, noticethat there the ratio is ∼ b and b before the actual CR would be greater than Anderson & Gordon C o un t s [ e - ] b b (a)DataFit to Semi-Ramps200002500030000350004000045000 C o un t s [ e - ] b b (b)2 4 6 8 10Sample Number020406080 | b (cid:13) b | / (cid:14) (c)With CRWithout CR Figure 5. y-int method. (a) Counts and two linear fits as itwould look when assuming that the CR is in sample 3. The solidline and points are the data, and the dotted lines are the linearfits to the semiramp before and after the assumed CR. The twoy-intercepts, b and b are highlighted with dots and labeled. (b)Assuming the CR is in sample 6 (which indeed it is). (c) Resultsafter stepping through the entire ramp assuming a CR is in eachsample. Plotted is the ratio of the absolute difference between they-intercepts to the uncertainty, which must be above the rejectionthreshold, r t , to be counted as a CR. The interval where the CRhit is shown by the vertical line. after the actual CR. RESULTSTo compare the success of each method at finding CRs,we look at both the fraction of CRs found and the falsedetection rate, which we define as the ratio of the num-ber of false detections to the number of possible falsedetections. The number of possible false detections iscalculated as the number of samples minus one, minusthe number of true CRs, times the number of trials.We look at simulated ramps with 5, 10, 20, 30, and40 samples, with CRs of 20 different magnitudes from 0to 250 e − (all CRs with magnitudes > e − are foundwith all methods), and with slopes of 0.0 e − /s, 0.7 e − /s,3.5 e − /s, 7.0 e − /s, 35.0 e − /s, and 70.0 e − /s, with one,two, or three CRs, and with the CR located in the be-ginning, middle, and end of the ramp (e.g. for a rampwith 40 samples, we inserted a CR in the 3 rd , 10 th , 20 th ,30 th , and the 35 th sample). We then simulate 10,000ramps for each combination, and apply each method toeach ramp to compare the results. For all trials, the re-jection threshold, r t , is chosen such that the rate of falsedetections is the same for all methods (a rate of 5% waschosen). In this way we can best compare how well theyfind the CRs.5.1. The Rejection Threshold, r t Each of the CR detection methods has a different cri-terion (see equations 10, 13, and 14), which, when com-pared with the rejection threshold, decides if a sampleis to be flagged as containing a CR or not. Therefore,the relationship between rejection threshold and the falsedetection rate is different for each method. To better un-derstand this dependence, for each CR rejection method F r a c t i o n o f F a l s e D e t e c t i o n s (b)0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5Rejection Threshold0.00.10.20.30.40.50.60.7 (c) m=0.0 e-/sm=0.7 e-/sm=3.5 e-/sm=7.0 e-/sm=35.0 e-/sm=70.0 e-/s Figure 6.
Fraction of false detections versus rejection threshold isplotted here for various input slopes. (a) method. (b) devfrom fit method. (c) y-int method. we used a range of rejection thresholds and looked at theresulting fraction of false detections. This was done onCR-free ramps.The first comparison is for the fraction of false detec-tions versus the rejection threshold for different inputslopes (Figure 6). For a given fraction of false detec-tions, the rejection threshold is independent of slope forthe method, while there is very little changefor the y-int method. However, for the dev from fit method, the rejection threshold needs to vary to keep thefraction of false detections constant across ramp slopes.The rejection threshold for different number of samplesin a ramp is shown in Figure 7. An input slope of 70.0 e − /s was used, and the rejection threshold was chosensuch that the fraction of false detections was 0.05. Theslope of these lines is 0.0011, 0.0065, and -0.0049 for the , dev from fit , and y-int methods respec-tively. Although this is a small change for the method, if a CR is found in a ramp then the rejectionthreshold was changed depending on the number of sam-ples in a semi-ramp. If the rejection threshold was notcalculated for a specific number of samples, we interpo-lated it. Figure 7 also shows how well the uncertaintiesin the method behave like a Gaussian, thusallowing the rejection threshold to be chosen easily.5.2. Slope and CR Detection
As we change the input slope of these ramps, we movefrom the read-noise-dominated regime to the photon-noise-dominated regime. In Figure 8 we show the twoextremes: read-noise-dominated regime with an inputslope of 0 e − /s (Figure 8 (b) and (d)), and photon-noise-dominated regime with an input slope of 70 e − /s (Fig-ure 8 (a) and (c)). Plotted is the fraction of CRs foundetecting Cosmic Rays in Infrared Data 7 R e j e c t i o n T h r e s h o l d Figure 7.
Rejection threshold that will result in a fraction of falsedetections of 0.05 is shown for various sample numbers. What wewould expect from a Gaussian distribution is also plotted to showhow similar it is to the method. A slope of 70.0 e − /s wasused for this figure. versus the CR magnitude and the fraction of false detec-tions versus the CR magnitude.In the photon-noise-dominated regime, photon-noisecan appear as a CR as it is correlated and thus can lead tofalse detections. From Figure 8 (a) and (c) we see that allof the methods give similar results. However, note howthe fraction of false detections for the methodis relatively constant compared with the other two meth-ods. This shows how well we are able to understand theuncertainties of this method. This, coupled with the factthat the method requires the least computa-tions, leads us to suggest that the method isthe optimal CR detection method in the photon-noise-dominated regime.In the read-noise-dominated regime we only have ran-dom noise, and therefore we are able to get a more ac-curate calculation of the y-intercept and the y-interceptuncertainties. Therefore, as is shown in Figure 8 (b)and (d), in the read-noise-dominated regime the y-int method gives the best results in that it is able to findfainter CRs than the other two methods.In Figure 8 (d), note that the rate of false detec-tions per CR strength is only constant for the method, while for the y-int method it is constant onlyfor high-energy CRs. Between a CR mag of 0 e − and 10 e − , we see a high fraction of false detections for the y-int method. This corresponds to the range where the CRsare not always found. The cause of this ’bump’ is notclear, and is evidence that we do not fully understandthe noise model associated with the Y-INT method.This could be caused by the fact that we treat eachsemiramp independently, which is inaccurate, since theuncertainties are correlated. Therefore, one solutioncould be to fit all semiramps simultaneously. This wouldincrease the uncertainty, which would mean that the re-jection threshold has been set too low (which would agreewith what we see in the plots). We leave that for futurework.Regarding Figure 8 (a) and (b), if you draw a verti-cal line at a given CR mag (we chose the CR mag wherethe fraction of CRs found by the method isclosest to 50%), you can see the difference between thefraction of CRs found by each method at that slope. We did this for slopes of 0.0 e − /s, 0.7 e − /s, 3.5 e − /s, 7.0 e − /s,35.0 e − /s, and 70.0 e − /s. The results are shown in Fig-ure 9. Each ramp had 40 samples, and again the rejec-tion threshold was chosen such that the fraction of falsedetections was 5%. For comparison, we subtracted thefraction of CRs found by the method. In thisfigure we can see that we move into the regime wherethe y-int method performs better than the and dev from fit methods at a slope of ∼ e − /s.To further illustrate how the results of each methodchange with input slope, we show the fraction of CRsfound as a function of CR magnitude for different slopesand for each of our methods in Figure 10. The results ofthe and dev from fit methods appear to besimilar, while the y-int method does better with lowerslopes (read-noise-dominated regime).5.3. Number of Samples and CR Detection
In order to determine how the number of samples ina ramp affects how well we detect CRs, we tested eachalgorithm on simulated ramps with 5, 10, 20, 30, and40 samples. Figure 11 shows the fraction of CRs foundas a function of CR magnitude for different numbers ofsamples for the , dev from fit , and y-int methods.When we use the method on weak CRs ( < − ) and strong CRs ( >
150 e − ), the fraction of CRs foundincreases with the number of samples (except for the casewhere there are five samples). For the dev from fit method, the fraction of CRs found increases with numberof samples for string CRs, but shows no difference forweak CRs. Finally, for the y-int method, the fractionof CRs found decreases with number of samples for weakCRs, but increases with number of samples for strongCRs. 5.4. CR Sample Number and CR Detection
There are two questions when it comes to CR samplenumber: the first is where do we find false detections inthe ramp, and the second is how well we find CRs ingiven positions in the ramp.To answer the first question, we created a histogramof the sample numbers of the false detections found onsimulated CR-free ramps for each method. These arepresented in Figure 12. The results show that both the and dev from fit methods are not biasedtoward any position in the ramp. However, the y-int method is biased toward the last samples in a ramp wherethere is less effect from the correlation between samples.To answer the second question, we simulated rampswith 40 samples and input slope of 70 e − /s, with oneCR each at sample (numbers 3, 10, 20, 30, or 35). Wethen applied each of the algorithms to these ramps andcompared the results, shown in Figure 13.The results of these plots all show that the is the only method that does not vary depending on thesample number of the CR, as can be expected. Both the dev from fit and y-int methods are weakly biased bysample number. This is also expected, since both requirefitting lines to data, and the fewer points in a line the lessaccurate the fit. Again, in the photon-noise-dominatedregime, the method is best.5.5. Multiple Cosmic Rays and CR Detection
Anderson & Gordon F r a c t i o n o f C R s F o un d m = 70.0 e-/s(a)0 50 100 150 200 250CR mag [e-]0.0460.0480.0500.0520.0540.0560.058 F r a c t i o n o f F a l s e D e t e c t i o n s m = 70.0 e-/s(c) m = 0.0 e-/s(b)0 10 20 30 40 50CR mag [e-]m = 0.0 e-/s(d) Figure 8.
CR detection rate and false detection rate as a function of CR magnitude for all three methods. In (a) and (c) the input slopewas 70 e − /s, in the photon-noise-dominated regime, where we recommend the method be used. In (b) and (d) the input slopewas 0 e − /s, in the read-noise-dominated regime, where the y-int method outperformed the other two. Each ramp has 40 samples. D i ff e r e n c e i n F r a c t i o n o f C R s F o un d Figure 9.
Difference between the fraction of CRs found by the dev from fit and y-int methods to the method. Fromthis plot we are able to see that the y-int method outperforms theother two with slopes less than ∼ e − /s. In a 1000 s integration, if we assume that 20% of thefield will be affected by CRs, then about 4% of the fieldcould be affected by two CRs, and 0.8% could be affectedby three. This is substantial enough that we need toaccount for this possibility.To test how well each method does at handling multipleCRs, we simulated 10,000 40-sample ramps for each of0-19 CRs. Each CR was of the same strength and inrandom but different samples. We show the results for 0-3 CRs in Figure 14.As can be seen, the method shows no differ-ence between one, two, or three CRs. We show only upto three CRs for simplicity, but we have found that the method does not start to show a difference inresults until about nine CRs and therefore should be theoptimal method in the photon-noise-dominated regime.Both the dev from fit and y-int methods, althoughshow more variance than the method. This iscaused by the fact that both methods require linear fitsthat will include two CRs for the y-int method (whenassuming the correct location of one of the CRs), and allthree CRs for the dev from fit method. Furthermore,for both methods the ramp will be segmented after thedetection of the first CR. DISCUSSIONThe three CR detection methods presented in this arti-cle were tested on simulated ramps, adjusted, corrected,and refined, in order to present optimum versions. Wenow compare our methods to determine which is bestsuited for various data.6.1.
Our results show that the method performsbest in the photon-noise-dominated regime. On top ofbeing straightforward and computationally simple, weshowed that the fraction of false detections versus rejec-etecting Cosmic Rays in Infrared Data 9 F r a c t i o n o f C R s F o un d (b) Dev. from Fit0 50 100 150 200 250CR mag [e-]0.00.20.40.60.81.0 (c) Y-Int. m=0.0 e-/sm=0.7 e-/sm=3.5 e-/sm=7.0 e-/sm=35.0 e-/sm=70.0 e-/s Figure 10.
CR detection rate as a function of CR magnitude forall methods. Each line is for a different input slope. There are 40samples in each ramp. F r a c t i o n o f C R s F o un d (b) Dev. from Fit0 50 100 150 200 250CR mag [e-]0.00.20.40.60.81.0 (c) Y-Int. Figure 11.
CR detection rate as a function of CR magnitude forall methods. Each line is for a different number of samples in aramp. The input slope is 70 e − /s for each ramp. O cc u rr e n c e Figure 12.
Histograms of the sample numbers of the false detec-tions for each method. These are from ramps without a simulatedCR. F r a c t i o n o f C R s F o un d (b) Dev. from Fit0 50 100 150 200 250CR mag [e-]0.00.20.40.60.81.0 (c) Y-Int. CR Sample
Figure 13.
CR detection rate as a function of CR magnitude forall methods. Each line is for a different CR sample number. Theinput slope is 70 e − /s, and there are 40 samples in each ramp.The method is the only one to not show a dependence onsample number. tion threshold is consistent regardless of slope and num-ber of samples. We also showed that the uncertaintiesof the method follow a Gaussian distribution.Any deviation from a Gaussian was initially thought tobe caused by the use of the median instead of averagein the rejection criterion. However we found that whileusing the average did produce a more Gaussian shape,it was still not a perfect Gaussian. It was demonstratedthat, as expected, the fraction of CRs found changedwith slope and number of samples (though barely). Fur-thermore, we showed that the false detections are evenlydistributed in all samples, and that the CR sample num-ber did not change the results. Finally, we found that0 Anderson & Gordon F r a c t i o n o f C R s F o un d (b) Dev. from Fit0 50 100 150 200 250CR mag [e-]0.00.20.40.60.81.0 (c) Y-Int. Figure 14.
CR detection rate as a function of CR magnitude.Each line is for a different number of CRs in a ramp. there is no noticeable difference in the performance ofthe method with multiple CRs, up to 9 CRson a ramp with 40 samples.Overall, the method is fast, consistent, andeasy to understand and calibrate (e.g., choose a rejectionthreshold), and it gives the best results in the photon-noise-dominated regime.6.2.
Y-Intercept Method
For the y-int method, there are two regimes:1. Photon-noise-dominated: In this regime the y-int method gives the same answer as the method. This is due to the correlated behaviorof the noise in the ramps (e.g., a 3 σ event due tothe noise of the photons that have arrived in thelast sample time will propagate through all sub-sequent samples, having the same effect as a CR).Calculating the linear fit of subsequent samples willnot negate the fact that the noise in the photon-dominated regime is set by the noise in a singlesample.2. Read noise-dominated: In this regime the y-int method is better than the method. Thesetwo methods will still give the same results on anynoise above the threshold, but anythingbelow that will only be detectable by the y-int method. With the noise independent from sampleto sample, a linear fit reduces the uncertainty inthe y-intercepts and weaker CRs are detectable.While it is unlikely that the uncertainty will be read-noise-dominated for most MIRI data (given the tele-scope and sky background), it will be the case for the
JWST NIRSpec , NIRCam , and the
Tunable Filter Im-ager ( TFI ). NIRSpec has a higher resolution than the
TFI and therefore a lower background. Regardless,the
TFI will likely be read-noise limited, given that ituses a tunable filter (narrowbands are less sky/telescope-dominated).The y-int method, however, does not quite match therobustness of the method. In theory, thesetwo methods should be the same in the photon-noise-dominated regime. However, the method gavemore consistent results for different slopes, number ofsamples, and CR location in the ramps and for multipleCRs. The difference is that the method has asimpler noise behavior, while we have not fully solved forthe noise behavior for the y-int method. Despite beingcareful when dealing with correlated and uncorrelatederrors when fitting a line to a ramp or semiramp, we didnot take this into consideration when taking the averageof the slopes of the semiramps. Instead we use a weightedaverage, where the weights are simply the uncertaintiesin the slope. A solution to this problem would be to fiteach semiramp and the CR simultaneously, as done byRobberto (2008) for uncorrelated errors only, modifyingthe method to account for correlated errors. Regardless,when working in the read-noise-dominated regime, the y-int method should be used.6.3.
Deviation from Fit Method
The dev from fit method did not perform as wellas the other two methods. The y-int method domi-nated in the read-noise regime, and in the photon-noise-dominated regime the dev from fit method, like the y-int method, was not as consistent as the method when it came to different slopes, number of sam-ples, CR location in the ramps, and multiple CRs. There-fore, to have a consistent fraction of false detections, wewould have to set a new rejection threshold based onthese variables. Furthermore, while not as complex atthe y-int method, the dev from fit method is stillmore complex than the method and thus com-putationally more expensive. Finally, the uncertaintiesare nowhere near Gaussian, thus making it difficult toset the rejection threshold. The dev from fit methodwas only slightly worse at detecting CRs than the method, but given the fact that it is not as robust,and requires more CPU time, the dev from fit methodis not recommended.The dev from fit method could be improved bychanging the weighting scheme so that instead of usingthe slope in the calculation of the photon-noise (sincethe slope will still include the CR in its calculation), wecould use the median of the two-point differences as wedo for the method. However, it still would notcompare with the method, as we cannot dobetter than two n r terms and one n p term by the natureof our problem. Furthermore, even if we were able to getthe same results for these two methods, the dev fromfit method still requires larger computational resourcesthan the method. CONCLUSIONSIn this article we discussed three CR detection methodsfor nondestructive read ramps, as well as their strengthsetecting Cosmic Rays in Infrared Data 11and weaknesses. We applied these methods to simulatedramps with single-sample groups. We showed that thecovariance matrix must include correlated errors in orderto improve the calculation of the slope and y-interceptof a ramp and their uncertainties. The y-int methodbenefits most from this use of the covariance matrix, andwould not be an improvement of the methodotherwise.The method was shown to be able to give thebest results in the photon-noise-dominated regime. Themethod’s uncertainties are quasi-Gaussian which simpli-fies the process of choosing a rejection threshold. More-over, it is fast and consistent. Its robustness, comparedwith the other two methods, resides in the fact that wesimply remove the two-point difference that includes theCR, instead of splitting the ramp into two semiramps,when searching for other CRs and calculating the slopeand the y-intercept.The dev from fit method results in a similar fractionof CRs found as the method for some ramps,but unlike the method these results changebased on slope, number of samples, and CR location inthe ramp and for multiple CRs. If we also consider thefact that it is computationally more expensive, we do notrecommend this method for use in any regime.The y-int method achieves the best results in theread-noise-dominated regime, and returns the same re-sults as the method in photon-noise-dominatedregime. In the photon-noise-dominated regime, the y-int method behaves like the method in thatit is effectively calculating an average slope, excludingdisturbances from any CRs (in the diff method this isdone by taking the median of the two-point differences).The average slope is then divided by the noise due toone sample time: one photon-noise and two read-noise.The y-int method is better in the read-noise-dominatedregime as the uncertainties on the line fit parameters di-minish in size as more points are fit. This is not the casein the photon-noise-dominated regime due to the corre-lated nature of the noise. The y-int method only hastwo drawbacks: the noise model is not fully understood,and it is a complex algorithm.In summary, if we take all results into consideration weare led to suggest that if only one method can be usedon data cubes, the y-int method should be it. If com-putational speed is an issue then the methodshould be used, especially where photon-noise dominates,and the y-int method should only be used when the es-timated slope is in the read-noise-dominated regime orwhere there is information, such as a nearby CR andcross-talk is suspected, to look for fainter effects. ACKNOWLEDGMENTSThis work was supported by the Space Telescope Sci-ence Institute, which is operated by AURA, Inc., underNASA Contract No. NAS5-26555.We thank Howard Bushouse for providing the informa-tion on the method used by NICMOS and WFC3, whichwas called the dev from fit method here.Appreciation goes to Harry Ferguson, David Grumm,Derck Massa, Mike Regan, and Massimo Robberto fortheir comments.2 Anderson & Gordonmethod here.Appreciation goes to Harry Ferguson, David Grumm,Derck Massa, Mike Regan, and Massimo Robberto fortheir comments.2 Anderson & Gordon