ramp_plot_4histo.pro (JWST Pipeline routines in IDL)

; docformat='rst'
;+
; This function 
;-    
   
;+
;  
;       
; :categories: 
;       JWST Data Processing - NIRCam Cryo tests
;       
; :examples:
;       IDL> RAMP_PLOT_4HISTO,frame,frame_4AVERAGES, MASK, file
;
; :Params: 
;   frame: in, required, type="fltarr(2048,2048)" 
;      frame to be plotted
;   frame_4averages: out, required, type="fltarr(5,4)"
;      the results of the gaussian fit
;   mask: out, required, type="fltarr(2048,2048)"   
;      the map of the bad pixels (>4sigma)   
;   fileinfo: in, optional, type="strarr(2)"
;     String array containing the PATH name [0], the directory name [1] and the filebasename [2] for postcript output
;      
; :Description:
;       gaussfit_results is a datacube with gaussfit results [5,2]::
;         intercept => [*,0]
;         slope     => [*,1]
;      Let us remember the 5 values returned by gaussfit,array,aa:: 
;      aa = coefficients of the Gaussian fit: Height, mean, sigma,...
;      
;         aa[0]= the height of the Gaussian::
;         aa[1]= the mean::
;         aa[2]= the standard deviation::
;         aa[3]= the half-width of the 95% conf. interval of the mean::
;         aa[4]= 1/(N-1)*total( (y-mean)/sigma)^2 ) = a measure of normality
;         
; :Author:
;       - Original from R. Arendt, adapted by M. Robberto, 25 Apr. 2013
;-
PRO RAMP_PLOT_4HISTO,frame,frame_4AVERAGES, MASK, fileinfo
st=systime(1)
print,'  begin RONMAP creation: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'

if N_ELEMENTS(fileinfo) NE 1 THEN BEGIN
   file_eps = fileinfo[0]+'/'+fileinfo[1]+'/'+fileinfo[2]+'_ron4histo.eps'
   ;
   ;check for the existence of the directory and if needed create it
   result = FILE_TEST (file_eps,/DIRECTORY) & if result EQ 0 THEN FILE_MKDIR, fileinfo[0]+'/'+fileinfo[1]
   ;open the .eps file
   openeps,file_eps
endif else begin 
  wset,0
  !p.multi = [0,2,2]
endelse

;====================================================================
;Masks outliers with especially high or low standard deviation
;Before moving forward with the analysis, we clean the data set. We look for outliers, with threshold set at 4 sigma

;Prepare the output 2x2 window 
;xx mask outliers with especially high or low stddev
  
;
;We now look at the four “quadrants”, or sectors of the array; their indexes have been previously labeled h1mask…h4mask. 
h=lindgen(2048,2048)
h1mask = h[0:511,*]
h2mask = rotate(h[512:1023,*],5)
h3mask = h[1024:512*3-1,*]     
h4mask = rotate(h[512*3:*,*],5) 
;downselect the active pixels; we don't care about the reference pixels here
h1mask=h1mask[4:*,4:2043]
h2mask=h2mask[*,4:2043]
h3mask=h3mask[*,4:2043]
h4mask=h4mask[4:*,4:2043]
;
;
;Plot the histogram of the standard deviations; the mean of the histogram is the readout noise, in electrons. 
;The histogauss function produces a plot, fits a Gaussian profile and returns the parameters of the fit.
; 
  sig_data=frame
;aspect_ratio=1.0
device,xsize=30,ysize=24
;DEVICE,/INCHES,XSIZE=10.0,YSIZE=8;,SCALE_FACTOR=0.5
  histogauss,sig_data[h1mask],aa1,charsize=1.5,_EXTRA={POS:[0.1, 0.55, 0.44, 0.93],xtitle:'sector 1',xticks:4}   ;intercept HISTOGAUSS is in the IDLASTRO library at http://idlastro.gsfc.nasa.gov
  histogauss,sig_data[h2mask],aa2,charsize=1.5,_EXTRA={POS:[0.6, 0.55, 0.94, 0.93],xtitle:'sector 2',xticks:4,noerase:1}   ;slope
  histogauss,sig_data[h3mask],aa3,charsize=1.5,_EXTRA={POS:[0.1, 0.07, 0.44, 0.45],xtitle:'sector 3',xticks:4,noerase:1}   ;slope
  histogauss,sig_data[h4mask],aa4,charsize=1.5,_EXTRA={POS:[0.6, 0.07, 0.94, 0.45],xtitle:'sector 4',xticks:4,noerase:1}   ;slope
  xyouts,0.05,0.95,fileinfo[2],charsize=1.8,/NORMAL
;
;Let us remember the values of aa
;      aa = coefficients of the Gaussian fit: Height, mean, sigma
;               aa[0]= the height of the Gaussian
;               aa[1]= the mean
;               aa[2]= the standard deviation
;               aa[3]= the half-width of the 95% conf. interval of the standard
;                     mean
;               aa[4]= 1/(N-1)*total( (y-mean)/sigma)^2 ) = a measure of 
;                       normality
;
;
;  histogauss,sig_data[h1mask],aa1   ; HISTOGAUSS is in the IDLASTRO library at http://idlastro.gsfc.nasa.gov
;  histogauss,sig_data[h2mask],aa2
;  histogauss,sig_data[h3mask],aa3
;  histogauss,sig_data[h4mask],aa4
if N_ELEMENTS(fileinfo) NE 1 THEN closeps

;
;Next we isolate the pixels that are within +/-4 sigma from the average readout noise: mask1=1 if the condition is satisfied, 
;0 otherwise. 
;Again, we extract from h1mask the indexes of the pixels with a discrepancy <4sigma (with respect to the average of the quadrant). 
;Mask 1 has values 1 or 0
  mask1 = abs(sig_data[h1mask]-aa1[1]) lt (4*aa1[2])
;Repeat for the other 3 sectors. There is here a “counter-reflection” around the y axis to have the pixels oriented normally; 
;this because h2mask and h4mask have been originally defined flipped with respect to the horizontal (fast readout) direction
  mask2 = rotate(abs(sig_data[h2mask]-aa2[1]) lt (4*aa2[2]),5)
  mask3 = abs(sig_data[h3mask]-aa3[1]) lt (4*aa3[2])
  mask4 = rotate(abs(sig_data[h4mask]-aa4[1]) lt (4*aa4[2]),5)
;idl returns 1 if the LT condition is satisfied, 0 otherwise. This procedure extracts from h1mask the values with  
;discrepancy <4sigma (with respect to the full array). Mask 1 has values 1 or 0
;




;
;THE RESULTS CAN BE EXTRACTED HERE
  frame_4AVERAGES=[[aa1],[aa2],[aa3],[aa4]]
;The variable mask contains now the indexes of all “good” pixels
  mask = [mask1,mask2,mask3,mask4]
 
 
 print,'    end RONMAP creation: ',systime(1)-st,' sec     ',memory(/high)/1e9,' GB'
END