;-
; docformat = 'rst'
;+
; This file contains the procedures to perform the reference pixel correction of NIRCam data
; according to the report by M. Robberto, 2014
;-
;+
; Performs referene pixel correction of NIRCam data
;
; :categories:
; JWST Data Processing - NIRCam Cryo tests
; :examples:
; IDL> corrector, ramp, rp_cleaned, /FILTER
; :Params:
; ramp: in, required, type=fltarr
; datacube of 2048,2048,N to be filtered
; rp_cleaned: out, required, type=float
; same datacube in input, after reference pixel correction
; /FILTER: in, optional, type=fltarr
; if set, the reference pixels are smoothed using the Kosarev and Pantos algorithm
; :uses:
; optimal_smooth_fft.pro
; :Author:
; - MR, February, 25, 2014
;-
pro corrector, ramp, rp_cleaned, FILTER = filter
;
CORRECTOR=ramp*0. ;same size of the input data, this will hold the corrector frame
;
;prepare the dataset of reference pixels
a = ramp ;do not overwrite the data, leave the ramp alone!
mean_counts = avg(a,2) ;average over all frames of the ramp
s = SIZE(ramp)
for i_fr=0,s[3]-1 do a[*,*,i_fr] -= mean_counts ;all pixels now bounce around their own average set at zero. Good only for ref.pixels.
;
;SECTOR1 and SECTOR4 provide the high frequency component.
;Start with SECTOR 1
sector1b = a[0:511,0:3,*] ;bottom reference pixels
sector1t = a[0:511,2044:2047,*] ;top reference pixels
sector1s = a[0:3,4:2043,*] ;side reference pixels ;side reference pixels
;
sector2b = reverse(a[512:1023,0:3,*],1)
sector2t = reverse(a[512:1023,2044:2047,*],1)
;
sector3b = a[1024:1535,0:3,*]
sector3t = a[1024:1535,2044:2047,*]
;
sector4b = reverse(a[1536:2047,0:3,*],1)
sector4t = reverse(a[1536:2047,2044:2047,*],1)
sector4s = reverse(a[2044:2047,4:2044,*],1) ;side reference pixels
;
;TIME is actually the same for all sectors (at least, the two outer sectors)
time = lindgen(512,2048)*10*1E-6 ;in seconds, assumes 10microsecond/pixel, not really critical but just to get the rates in good units
timeb = time[0:511,0:3] ;readout time of the bottom reference pixels
timet = time[0:511,2044:2047] ;readout time of the top reference pixels
times = time[0:3,4:2043] ;readout time of the side reference pixels
;
; LOOP ON THE FRAMES
for i_fr=0,s[3]-1 do begin
;
;the linear fit is performed only using bottom/top refpix., to be consistent with the central sectors.
;1) find the linear fit to the top/bottom ref.pix.
coeffs1=linfit([timeb[*,*],timet[*,*]],[sector1b[*,*,i_fr],sector1t[*,*,i_fr]])
coeffs4=linfit([timeb[*,*],timet[*,*]],[sector4b[*,*,i_fr],sector4t[*,*,i_fr]])
;
;extract the high frequency signal
;use ALL reference pixels, bottom and top, to construct a 2048 vector
rp1 = [avg(sector1b[*,*,i_fr],0),avg(sector1s[*,*,i_fr],0),avg(sector1t[*,*,i_fr],0)]
time_rp1 = [avg(timeb[*,*],0),avg(times[*,*],0),avg(timet[*,*],0)]
;
;the high frequency term is calculated taking the vector of 2048 row averaged RefPixels
;(either 512 or 4, so big difference here...) and subtracting the slote derived by the use of the bottom/top
;reference pixels only. So, basically, adding back to this HIGHFREQ vector the slope I return to the
;averaged RefPix values, working ecumenically for all sectors as long as I get the stopes in the
;same way, from the bottom/top refpix.
highfreq_1=rp1-coeffs1[0]-coeffs1[1]*time_rp1 ;all reference pixels are corrected for the linear fit
;
;repeat for SECTOR4
;find the fit line
rp4 = [avg(sector4b[*,*,i_fr],0),avg(sector4s[*,*,i_fr],0),avg(sector4t[*,*,i_fr],0)]
highfreq_4=rp4-coeffs4[0]-coeffs4[1]*time_rp1
HF = (highfreq_1+highfreq_4)/2.
;
;;HERE WE CAN FILTER THE VECTOR...
IF keyword_set(SILENT) THEN BEGIN
optimal_smooth_fft, HF, 5E-3, Smoothed_Data
HF=smoothed_Data
ENDIF
;
;now move to the central sectors
;SECTOR2
coeffs2=linfit([timeb[*,*],timet[*,*]],[sector2b[*,*,i_fr],sector2t[*,*,i_fr]])
;
;SECTOR3
coeffs3=linfit([timeb[*,*],timet[*,*]],[sector3b[*,*,i_fr],sector3t[*,*,i_fr]])
;
print,i_fr,coeffs1,coeffs2,coeffs3,coeffs4
;
;reconstruction
corrector[0:511,*,i_fr]=(coeffs1[0]+coeffs1[1]*time_rp1+HF)##replicate(1,512)
corrector[512:1023,*,i_fr]=(coeffs2[0]+coeffs2[1]*time_rp1+HF)##replicate(1,512)
corrector[1024:1535,*,i_fr]=(coeffs3[0]+coeffs3[1]*time_rp1+HF)##replicate(1,512)
corrector[1536:2047,*,i_fr]=(coeffs4[0]+coeffs4[1]*time_rp1+HF)##replicate(1,512)
endfor
;
;and the Residual Ref.Pix.cleaned frame is
rp_cleaned=ramp-corrector
;
;LAST STEP: PEDESTAL CORRECTION
;==============================
RP_mask = rp_cleaned
RP_mask[4:2043,4:2043,*] = !VALUES.F_NAN
RP_mean = mean(RP_MASK,/NAN)
RP_meanS = fltarr(4,s[3])
RP_meanS[0,*] = mean(RP_MASK[0:511,*,*],/NAN)
RP_meanS[1,*] = mean(RP_MASK[512:1023,*,*],/NAN)
RP_meanS[2,*] = mean(RP_MASK[1024:1535,*,*],/NAN)
RP_meanS[3,*] = mean(RP_MASK[1536:2043,*,*],/NAN)
;
for i_fr=0,s[3]-1 do begin & $
rp_cleaned[0:511,*,i_fr] += (RP_mean-RP_meanS[0,i_fr]) & $
rp_cleaned[512:1023,*,i_fr] += (RP_mean-RP_meanS[1,i_fr]) & $
rp_cleaned[1024:1535,*,i_fr] += (RP_mean-RP_meanS[2,i_fr]) & $
rp_cleaned[1536:*,*,i_fr] += (RP_mean-RP_meanS[3,i_fr]) & $
endfor
end