;pro PlotSlope ; ;a routine to plot the slopes of pixels for diagnostic or presentation ;purposes pro PlotSlope, Col, Row, Raw, Dark, TypeString Slope = Raw - Dark plot, Slope, yrange=[min(Slope), max(Slope)],$ title = 'Row : '+Strtrim(Row,2) + ', Col : '+Strtrim(Col,2) stop end ;pro SlopeFit ; ;PURPOSE: ; ; To fit a slope to each of the pixels in a datacube. The slope of ; the pixel value should represent the signal seen by that pixel, ; whether it is darkcurrent or light. ; ; This routine should properly handle saturated pixels in a multi-read ; exposure with the exception of pixels that got saturated on the first read. ;NOTES: ; This should accomodate both raw images and median dark frames if the ; filenames are set correctly. That is, if the file name contains the ; word 'Bias' or 'Dark', there will be no attempt to subtract a dark ; frame. ; ;KEYWORDS: ; TVFLAG: 0 - Default, Don't plot anything ; 1 - Plot to the X window ; 2 - Plot to postscript ; 3 - Plot to a .png ; ; TVPPT: Tv Plot Pixel Type - Designates the type of pixel to be plotted ; 1 - Plot ; 2 - Plot the first and second read saturated pixels ; 3 - Plot the cosmic rays that occur somewhere in the ramp ; 4 - Plot the first read cosmics or potentially hot pixels ;CALLING SEQUENCE ; ASIC Images--- ; slopefit_1rg,objectname='M1',date='07Nov17',thisfilter='G', /subtractdark,/userefpix, nreads=250 ; LEACH Images---- ; Individual Image ; slopefit_1rg,'E:\LEACH\2007Nov18\Orion_Trapezium\Orion_Trapezium_Gfilter_15_1_0_0_0_VSUB15_0.fits',/subtractdark,date='2007Nov18',electronics=1 ; Batch Images:; ; slopefit_1rg,objectname='Orion_Trapezium',/subtractdark,date='2007Nov18',electronics=1, objecttext='Gfilter_15_1_0_', /userefpix, thisfilter='G' ; slopefit_1rg,objectname='Orion_Trapezium',/subtractdark,date='2007Nov18',electronics=1, objecttext='Ifilter_15_1_0_', /userefpix, thisfilter='I' ; ;REFRENCES: ; Based on code used in CalnicA for the NICMOS Data Reduction ;**************************************************************************** pro Slopefit_1RG, RawFileName, RootDir = RootDir, Date = Date, $ ReducedDir = ReducedDir, SubtractDark = SubtractDark, $ ObjectName = ObjectName, ObjectText=ObjectText, $ NReads = NReads, ThisNDrops = ThisNDrops, $ TvFlag = TvFlag, TvPPT = TvPPT, UpdateMask=UpdateMask, $ ThisFilter = ThisFilter, SkipFileNum=SkipFileNum, $ UseRefPix = UseRefPix, Offset = Offset, UseBadMask = UseBadMask, $ Electronics=Electronics If N_Elements(UpdateMask) eq 0 then UpdateMask = 0 If N_Elements(SkipFileNum) eq 0 then SkipFileNum = 0 If N_Elements(TvPPT) eq 0 then TvPPT = 0 If N_Elements(Electronics) eq 0 then Electronics = 0 If N_Elements(UseRefPix) eq 0 then UseRefPix = 0 Common KeyParams, KeyStr ;Include all of the keywords from the KeywordStruct file @KeywordStruct_1RG.pro @PlotSettings_1RG.pro ;************************************************************************* ;Move to the directory If stregex(RootDir, 'Reduced', /boolean) then begin cd, RootDir Endif Else begin If Electronics eq 0 then begin cd, RootDir + 'ASIC'+KeyStr.PathDelim+Date EndIf else begin cd, RootDir + 'LEACH'+KeyStr.PathDelim+Date EndElse EndElse ;******BATCH MODE********************************************************* ;If filename is blank, search for matching files in the directory If N_Elements(RawFileName) eq 0 then begin If KeyStr.Electronics eq 0 then begin RawFileName = File_Search(KeyStr.ObjectName+'*'+$ strtrim(KeyStr.NReads, 2) + '_Reads'+'*.fits') RawFileName = Return_File_List_1rg(RawFileName, ThisFilter=ThisFilter, $ ThisNDrops=ThisNDrops) EndIf Else Begin cd, KeyStr.ObjectName RawFileName = File_Search(KeyStr.ObjectName+'*'+$ KeyStr.ObjectText+'*.fits') EndElse RawFileName = RawFileName(SkipFileNum:N_Elements(RawFileName)-1) print, 'About to calculate slopes for: ' print, RawFileName KeyStr.ObjectDir = KeyStr.ReducedDir+KeyStr.ObjectName+KeyStr.PathDelim File_MkDir, KeyStr.ObjectDir EndIf ;**************************************************BEGIN LOOP OVER FILES ;Batch mode for given filter or single image For FileNum = 0, N_Elements(RawFileName)-1 do begin print, 'Working on File :' + RawFileName(FileNum) ;Get the current file that's being operated on and some keys KeyStr.RawName = RawFileName(FileNum) ;Get the Parameters necessary to fit this file If Electronics eq 0 then begin Fits_Read, KeyStr.RawName, NoData, RawHeader, /header_only FileParam= Return_File_Params_1RG(RawFileName(FileNum)) SlopeStr = Return_Slope_Params_1RG(RawHeader) EndIf Else Begin FileParam= Return_LEACH_File_Params_1RG(RawFileName(FileNum)) SlopeStr = Return_LEACH_Slope_Params_1RG(KeyStr.RawHeader) EndElse ;Check if File is a dark itself, if it is, don't subtract and If stregex(KeyStr.RawName, 'Median', /boolean) and $ (stregex(KeyStr.RawName, 'Dark', /boolean) or $ stregex(KeyStr.RawName, 'Bias', /boolean)) then begin KeyStr.UseBadMask =0 KeyStr.SlopeSigThresh =6 EndIf ;Make a temporary mask that can be chopped to save time BadInd = 0L NumNonLin = 0L NReads = KeyStr.TotalReads TimePoints = SlopeStr.ReadTimes(1:NReads-1) LowerLinLim = KeyStr.LowerLinLim UpperLinLim = KeyStr.UpperLinLim FullLinLim = UpperLinLim - LowerLinLim ;Make arrays for the Slopes and what have you Slopes = Fltarr(KeyStr.Naxis1, KeyStr.Naxis2) Intercepts = Fltarr(KeyStr.Naxis1, KeyStr.Naxis2) Errors = Fltarr(KeyStr.Naxis1, KeyStr.Naxis2) Cosmics = Fltarr(KeyStr.Naxis1, KeyStr.Naxis2) SatMask = BytArr(KeyStr.Naxis1, KeyStr.Naxis2) ;Attempt to remove fixed noise using channel 0 If KeyStr.UseRefPix then begin SmooRefCols = Return_RefPix_1RG(KeyStr.RawName,0) SmooRefColL = (SmooRefColR = Fltarr(KeyStr.Naxis2,KeyStr.Naxis3)) SmooRefColL(*,*) = Float(SmooRefCols(0,*,*)) SmooRefColR(*,*) = Float(SmooRefCols(1,*,*)) EndIf ;Take in Frames one half at a time to avoid memory crash For FrHalf = 0, 1 do begin print, 'Computing Slopes for frame half : '+strtrim(FrHalf,2) stop ;Read in Half of the raw/median image Fits_Read_DataCube, KeyStr.RawName , RawHalf , RawHeader, $ XStart= 0, XStop = KeyStr.Naxis1-1, $ YStart= FrHalf*KeyStr.Naxis2/2, $ YStop = (FrHalf+1)*KeyStr.Naxis2/2-1, $ Offset = Offset stop ;Read in Half of the dark image if dark subtraction is requested If KeyStr.SubtractDark then $ Fits_Read_DataCube, KeyStr.DarkName, DarkHalf, DarkHeader, $ XStart= 0, XStop = KeyStr.Naxis1-1, $ YStart= FrHalf*KeyStr.Naxis2/2, $ YStop = (FrHalf+1)*KeyStr.Naxis2/2-1 $ Else DarkHalf = LonArr(KeyStr.Naxis1, KeyStr.Naxis2, KeyStr.Naxis3) ;Loop over the pixels and slope fit them individually and iteratively, ;rejecting cosmic rays For Row = Double(FrHalf*KeyStr.Naxis2/2), $ Double((FrHalf+1)*KeyStr.Naxis2/2 - 1) do begin ;Get the index for the half frame RowH = Row-FrHalf*KeyStr.Naxis2/2 print, Row ;Remove the common mode from each row in each frame For Col = 0D, KeyStr.Naxis1-1 do begin PixelInd = Col+KeyStr.Naxis1*Row If KeyStr.SubtractDark then begin DarkRamp = DarkHalf(Col,RowH,1:KeyStr.TotalReads-1) EndIf Else begin DarkRamp = LonArr(KeyStr.TotalReads) EndElse ;Take the Raw Ramp in to check for Saturation without any ;subtration taking place RawRamp = RawHalf(Col,RowH,1:NReads-1) RawBias = RawHalf(Col,RowH,0) If Not KeyStr.DarkSlope then begin ;Find values in the linear regim GoodMask = where(RawRamp gt LowerLinLim and $ RawRamp lt UpperLinLim, NumGood) DarkFirstRead = DarkHalf(Col,RowH,1) DarkBias = DarkHalf(Col,RowH,0) RawFirstRead = RawHalf(Col,RowH,1) ;Check for saturated pixels and add them to the saturation mask SatVals = where( RawRamp gt UpperLinLim, NumSat) If SatVals(0) ne -1 then SatMask(Col,Row) = NumSat ;Subtract the Reference Ramp after the limits of linearity ;have been checked If KeyStr.UseRefPix eq 1 then begin RawRamp = RawRamp - SmooRefCols(Col/(KeyStr.Naxis1/2), RowH, $ 1:KeyStr.TotalReads-1) RawBias = RawBias - SmooRefCols(Col/(KeyStr.Naxis1/2), RowH, 0) EndIf EndIf Else Begin GoodMask = where( (RawRamp +RawBias) gt LowerLinLim and $ (RawRamp +RawBias) lt UpperLinLim, NumGood) DarkFirstRead = 0 DarkBias = KeyStr.LowerADCThreshold RawFirstRead = RawHalf(Col,RowH,1)+RawHalf(Col,RowH,0) EndElse ;Perform a line fit to the points in the linear region If NumGood le 1 then begin NumNonLin +=1 ;If pixel is saturated, check neigbhors to see if it's at the ; center of a star. If so, use read 0 as an estimate of flux ;Make sure that borders of array are not a problem ColLeft = 1 & ColRight = 1 & RowDown = 1 & RowUp = 1 If Col eq 0 or Col eq KeyStr.Ch0W then ColLeft = 0 If Col eq KeyStr.Naxis1 - 1 then ColRight = 0 If RowH eq KeyStr.Ref/2 or RowH eq 0 then RowDown = 0 If RowH eq KeyStr.Naxis2/2-1 then RowUp = 0 ;Get the 3x3 (or smaller) region around RegCheck = RawHalf(Col-ColLeft:Col+ColRight,$ RowH-RowDown:RowH+RowUp,0) ;Get the counts of the flagged pixel CenFlux = RegCheck(ColLeft,RowDown) ;Set it to zero and look at the avg. value of neighbors RegCheck(ColLeft,RowDown) = 0 SurrPix = where(RegCheck ne 0, NumSurrPix) If NumSurrPix ne 0 then begin djs_iterstat, RegCheck(where(RegCheck ne 0)), $ mean = SurMean, sigma = SurSigma Endif Else begin SurMean = 0 & SurSigma = 0 EndElse ;If it's an outlier, flag it as a cosmic hit or a bad pixel ;that hasn't been discovered yet If CenFlux gt SurMean+2.5*SurSigma and $ Not KeyStr.DarkSlope then begin Slopes(Col,Row) = 0 print, 'Hot Pixel or Cosmic at : ', Col, Row If KeyStr.TvFlag and KeyStr.TvPPT eq 4 then $ PlotSlope, Col, Row, RawHalf(Col,Row,*), DarkHalf(Col,Row,*), $ 'Hot Pixel' EndIf else begin ;If the first read is still in the linear regime and it's ;less than half the way there, use the single cds If NumGood eq 1 and $ RawBias lt (float(FullLinLim)/2+LowerLinLim) $ then begin Slopes(Col,Row)= $ ((RawFirstRead-RawBias)-$ DarkFirstRead)/$ (SlopeStr.GroupTimes(1)-SlopeStr.GroupTimes(0)) EndIf Else begin ;If the first read is already saturated, subract ;the median dark bias and divide by a frame time Slopes(Col,Row)=(RawBias-DarkBias)/$ SlopeStr.FrameTime EndElse print, 'SatPix approx slope = ', Slopes(Col,Row) EndElse EndIf else begin LineFit_1RG, TimePoints, RawRamp, m, b, $ sigm, sigb, Error, CosmicEDep, $ Y2=DarkRamp, Mask=GoodMask, Col=Col, Row=Row, $ ThresholdFactor = KeyStr.SlopeSigThresh, $ TvFlag=KeyStr.TvFlag, TvPPT=KeyStr.TvPPT Slopes(Col,Row) = m Intercepts(Col,Row) = b Errors(Col,Row) = Error Cosmics(Col,Row) = CosmicEDep EndElse EndFor EndFor EndFor ;=========================================================================== ;*************************************************************************** ;SLOPEFIT OUTPUT FITS FILE ************************************************* ;0 = Slope fit ;2 = Error in Fit ;1 = Y Intercept ;3 = Cosmic Rays ;Form a header with some of the necessary keywords SlopeArray=[[[Slopes]],[[Intercepts]],[[Errors]],[[Cosmics]]] MkHdr, SlopeRawHeader, SlopeArray SlopeRawHeader = Update_Header_1RG(SlopeRawHeader) FxAddPar, SlopeRawHeader, 'FRAME1', 'Slope', 'm in y=mx+b' FxAddPar, SlopeRawHeader, 'FRAME2', 'Intercept', 'b in y=mx+b' FxAddPar, SlopeRawHeader, 'FRAME3', 'Error', $ 'Error in fit sqrt(<(data - y=mx+b)^2>)' FxAddPar, SlopeRawHeader, 'MINTHRESH', KeyStr.LowerLinLim, $ 'Minimum limit for linear slope calculation' FxAddPar, SlopeRawHeader, 'MAXTHRESH', KeyStr.UpperLinLim, $ 'Maximum limit for linear slope calculation' FxAddPar, SlopeRawHeader, 'NOISEREM', KeyStr.Rem60HzStr, $ 'Mode of noise removal for raw image' FxAddPar, SlopeRawHeader, 'MEDDARK', KeyStr.DarkKey, 'Name of Dark used' FxAddPar, SlopeRawHeader, 'NOISEREM', KeyStr.Rem60HzStr, $ 'Mode of fixed noise removal' FxAddPar, SlopeRawHeader, 'PROFILE', 'slopefit_1rg.pro', $ 'IDL File used to create slope' ;Form the filename depending on whether is was a dark slope or a lit slope ;UseRefPix should already be included in string If KeyStr.DarkSlope ne 1 then begin If KeyStr.SubtractDark eq 1 and $ Stregex(KeyStr.RawKey, 'NoDarkSubtracted', /boolean) then begin KeyStr.RawKey = strmid(KeyStr.RawKey, 0 , $ strpos(KeyStr.RawKey, 'NoDarkSubtracted')) KeyStr.RawKey = KeyStr.RawKey+'_' EndIf SlopeFileName = KeyStr.ObjectDir + KeyStr.RawKey+KeyStr.ThisFilter+$ '_Slope.fits' EndIf else begin ;Get rid of the 'NoDarkSubtracted' string If KeyStr.SubtractDark eq 1 and $ Stregex(KeyStr.RawKey, 'NoDarkSubtracted', /boolean) then begin KeyStr.RawKey = strmid(KeyStr.RawKey, 0 , $ strpos(KeyStr.RawKey, 'NoDarkSubtracted')) KeyStr.RawKey = KeyStr.RawKey+'_' EndIf SlopeFileName = KeyStr.DarkSlopeDir + KeyStr.Date+'_'+$ KeyStr.RawKey+KeyStr.ThisFilter+$ '_Slope.fits' EndElse Fits_Write, SlopeFileName, SlopeArray, SlopeRawHeader print, 'Slopes written to : ' , SlopeFileName ;***********************************************FLAT FIELD SLOPE ;1RG Files use the slopefitted flats with a certain amount of reads ;Divide by the flat if requested If not KeyStr.DarkSlope then begin If KeyStr.FlatName ne '' and $ Not Stregex(KeyStr.ObjectName, 'Flat', /boolean) and $ Not Stregex(KeyStr.ObjectName, 'Dark', /boolean) then begin Fits_Read, KeyStr.FlatName, FlatIm, FlatHeader ;Normalize the Flat FlatIm = Double(FlatIm(*,*,0))/Max(FlatIm(*,*,0)) FlatFieldSlope = Double(SlopeArray(*,*,0))/Double(FlatIm) FxAddPar, SlopeRawHeader, 'FLATFILE', KeyStr.FlatKey, $ 'Flat field used for image' FlatSlopeFileName = KeyStr.ObjectDir+KeyStr.RawKey+KeyStr.ThisFilter+$ KeyStr.Rem60HzStr+'_SlopeFlatFielded.fits' Fits_Write, FlatSlopeFileName, $ FlatFieldSlope, SlopeRawHeader Fits_Open, FlatSlopeFileName, SlopeFCB,/append Fits_Write, SlopeFCB, SatMask, extname = 'SatMask' Fits_Close, SlopeFCB EndIf Endif ;**************************************************************** EndFor stop End