;pro Return_Pixel_History ; ;this routine will get all of the values of a given pixel from the nine dithers ;Get the value of the saturated pixels for all of the dithers ;Read a pixel from file at a time (seems to be faster than reading ;entire datacube for less than 10,000 pixels or so pro Return_Pixel_History, X, Y, Times, PixHstry, PixMDarkHstry, $ ReadTimes, Slopes, Intercepts Common KeyParams, KeyStr TReads = KeyStr.NReads NumRawFiles = KeyStr.NumDithers RawFileNames = KeyStr.DitherFiles(0:NumRawFiles-1) ;Set up arrays to hold the various quants PixHstry = UIntArr(NumRawFiles*TReads) PixMDarkHstry = LonArr(NumRawFiles*TReads) ReadTimes = FltArr(NumRawFiles*TReads) Slopes = Fltarr(NumRawFiles) Intercepts = Fltarr(NumRawFiles) ;Get the dark for subtraction and fitting purposes Fits_Read_DataCube, KeyStr.DarkName, Dark, $ XStart = X, XStop = X, $ YStart = Y, YStop = Y, $ ZStart = 0, ZStop = TReads-1 ;For this pixel, go through the files For FileNum = 0, NumRawFiles - 1 do begin Fits_Read_DataCube, RawFileNames(FileNum), Im, $ Offset = KeyStr.Offset, $ XStart = X, XStop = X, $ YStart = Y, YStop = Y, $ ZStart = 0, ZStop = TReads-1 ReadTimes(FileNum*TReads:$ (FileNum+1)*TReads-1)=$ Times(FileNum)+Indgen(TReads)*KeyStr.ExpTime-$ Times(0) PixHstry(FileNum*TReads:$ (FileNum+1)*TReads-1)=$ Im(0,0,*) PixMDarkHstry(FileNum*TReads+1:$ (FileNum+1)*TReads-1)=$ Long(Im(1:TReads-1))-$ Long(Dark(1:TReads-1)) ;Do a line fit GoodMask = where(Im(1:TReads-1) gt KeyStr.LowerLinLim and $ Im(1:TReads-1) lt KeyStr.UpperLinLim, NumGood) If NumGood gt 2 then begin LineFit, (Indgen(TReads-1)+1)*KeyStr.ExpTime, $ Im(1:TReads-1), Dark(1:TReads-1),$ GoodMask, m, b, sigm, sigb, Error, $ CosmicEDep, X, Y, KeyStr.SlopeSigThresh, 0, 0 Slopes(FileNum) = m Intercepts(FileNum) = b EndIf EndFor End ;pro Examine_Persistent_Pixels ; ;PURPOSE: ; ; To take in a mask of pixels that were saturated in an image and ; plot the behavior of those pixels in the images immediately following. ; In the case of the data taken at KPNO, we are looking at a sequence of ; dithers each of which are separated in time by less than a minute. ; ;INPUTS: ; SlopeFileName - The name of the slopefit that contains the mask ; of saturated pixels. In most cases it should ; be the 0th dither ; ;KEYWORDS: ; ; MODE - ; 0 - Plot 1 set of slopes for each radius for one star ; 1 - Plot the average of the slopes vs. radius for each star ; 2 - Scatter plot the values for a given read number vs. dither ; ;CALLING SEQUENCE: ; Examine_Persistent_Pixels, '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr29/M16/M16_Dither0_0ArcSecondsEast_0ArcSecondsNorth_H4RG_SIPIN_15_Reads_Apr30_2007_03_48_59I_SlopeFlatFielded.fits', Offset=1440 ; ; Examine_Persistent_Pixels, '/nfs/slac/g/ki/ki09/lances/Reduced/07Apr29/M16/M16_Dither0_0ArcSecondsEast_0ArcSecondsNorth_H4RG_SIPIN_20_Reads_Apr30_2007_04_10_57G_SlopeFlatFielded.fits', Offset=1440, mode=1 ;**************************************************************************** pro Examine_Persistent_Pixels, SlopeFileName, $ RootDir = RootDir, Date = Date, $ ReducedDir = ReducedDir, ObjectName = ObjectName, $ SubtractDark = SubtractDark, FilterStr = FilterStr, $ TvFlag = TvFlag, Rem60Hz = Rem60Hz, $ Offset = Offset, Mode = Mode If N_Elements(Mode) eq 0 then Mode = 0 Common SlopeParams, SlopeStr Common KeyParams, KeyStr Common ColorStr, Color, ColorNum, ColorTri, ColorNumTri ;Include all of the keywords from the KeywordStruct file @KeywordStruct.pro @PlotSettings.pro ;Get the saturated mask from the slope file FileParams = Return_File_Params(SlopeFileName, /slope) Fits_Read, SlopeFileName, SatMask, exten_no = 1 ;Get the Files that are taken immediately after the first dither RawFileNames = File_Search(KeyStr.RootDir+KeyStr.PathDelim+$ KeyStr.Date+KeyStr.PathDelim+KeyStr.ObjectName+'*'+$ strtrim(KeyStr.NReads, 2) + '_Reads'+'*.fits') RawFileNames = Return_File_List(RawFileNames, ThisFilter=KeyStr.ThisFilter, $ ThisNDrops=ThisNDrops, Times=Times) NumRawFiles = N_Elements(RawFileNames) KeyStr.NumDithers = NumRawFiles KeyStr.DitherFiles(0:NumRawFiles-1) = RawFileNames print, RawFileNames ;Determine the name of the output file If KeyStr.TvFlag eq 2 then $ device, filename = KeyStr.PlotDir+KeyStr.ObjectName+'_RadialPersistence.ps' ;prepare some plotting indices TReads = KeyStr.TotalReads TReads = KeyStr.NReads PlotInds = Indgen(NumRawFiles*TReads) PlotInds = PlotInds(where(PlotInds mod TReads ne 0)) PlotInds = PlotInds(where(PlotInds gt TReads)) Masks = Return_Pixel_Masks() ;Find the regions that were occupied by bright stars as opposed to ;very high E cosmic rays or messed up pixels FitRad = 6. HMin = 2. ; minimum counts necessary to see a peak for find.pro FWHM = 6. ; FWHM cutoff for find.pro RoundLim = [-2.,2.] ; roundness paramtetrs for find.pro SharpLim = [0.2,1.] ; sharpness parameters for find.pro PhotAper = [2.] ; aperture radius used by aper.pro RdNoise = [3.] ; readnoise used by psf.pro PhPerADU = [1.] ; number of photons per adu PSFRad = 20. ; radius of the circular area in which psf will be def FitRad = 6. ;IRAF like routines to find sky, stars, etc. group will remove doubly counted Sky, SatMask, SkyMode, SkySig Find, SatMask, XStar, YStar, Flux, Sharp, Roundness, HMin, $ FWHM, RoundLim, SharpLim,/silent Aper, SatMask, XStar, YStar, Mag, ErrAp, Sky, SkyErr, 1., PhotAper, -1, $ [-1., TReads+1], setskyval = SkyMode, /silent Group, XStar, YStar, 20, GroupNum ;These are the approximate centroids of the stars that will be used XStar = XStar(Uniq(GroupNum)) YStar = YStar(Uniq(GroupNum)) NumStars = N_Elements(XStar) ;************************************************************************** ;The stars have been found; Now use them in one of the following modes Case Mode of ;************************************************Individual Slopes 0: begin XSigs = Fltarr(NumStars) & YSigs = Fltarr(NumStars) StarNum = 0 XStart = Round(XStar(StarNum)) - 25 > 0 XStop = Round(XStar(StarNum)) + 25 < KeyStr.Naxis1 YStart = Round(YStar(StarNum)) - 25 > 0 YStop = Round(YStar(StarNum)) + 25 < KeyStr.Naxis1 NewStar: Centroid = Find_Centroid_Sat(SatMask, XStart, XStop, YStart, YStop, $ IsStar, SatRad, TvFlag = KeyStr.TvFlag) XSigs(StarNum) = Centroid(2) & YSigs(StarNum) = Centroid(3) ;Now the centroid is confirmed and it's a star, look at the radii SatRegSize = 101 & RadMax = sqrt(2*(SatRegSize/2)^2) RadMax = 50 ;Check if we'll stay in bounds If Centroid(0) + SatRegSize/2 gt KeyStr.Naxis1 or $ Centroid(0) - SatRegSize/2 lt 0 or $ Centroid(1) + SatRegSize/2 gt KeyStr.Naxis2 or $ Centroid(1) - SatRegSize/2 lt 0 then goto, NewStar RelXCoords = Indgen(SatRegSize,SatRegSize) mod SatRegSize RelXCoords = RelXCoords - Mean(RelXCoords) RelYCoords = reverse(Indgen(SatRegSize,SatRegSize),2) / SatRegSize RelYCoords = RelYCoords - Mean(RelYCoords) RsqCoords = RelXCoords^2+RelYCoords^2 ;Plot one set of slopes for each radius For Rad = 0, RadMax - 1 do begin print, 'At a radius of : ' +Strtrim(Rad,2) RelPixInRad = where(sqrt(RsqCoords) ge Rad $ and sqrt(RsqCoords) lt Rad+1) AbsPixInds = RelXCoords(RelPixInRad) + Centroid(0) + $ (RelYCoords(RelPixInRad) + Centroid(1))* KeyStr.Naxis1 FirstPixel = 1 For PixNum = 0, N_Elements(AbsPixInds) - 1 do begin XInd = long(AbsPixInds(PixNum)) mod KeyStr.Naxis1 YInd = long(AbsPixInds(PixNum)) / KeyStr.Naxis2 If Masks.BadMask(XInd, YInd) eq 1 then break Return_Pixel_History, XInd, YInd, Times, PixHstry, $ PixMDarkHstry, ReadTimes, Slopes, Intercepts ;Now do the plotting IF KeyStr.TvFlag gt 0 and FirstPixel eq 1 then begin !p.multi=[0,1,2] Plot, ReadTimes, PixHstry, $ Title = 'Pixel History for ('+$ Strtrim(XInd,2)+','+Strtrim(YInd,2)+')', $ xtitle = 'Time (seconds)', ytitle = 'Raw Counts', $ background = fsc_color('white'), color=fsc_color('black') Xyouts, 0.7,0.9, 'Radius from Centroid : ' +Strtrim(Rad,2), $ color=fsc_color('black'), /normal Plot, ReadTimes, PixHstry, /nodata, $ yrange=[Min(PixMDarkHstry(PlotInds)),$ Max(PixMDarkHstry(PlotInds))],$ xtitle = 'Time (seconds)', ytitle = 'Raw-Dark Counts', $ color=fsc_color('black') For FileNum = 1, NumRawFiles - 1 do begin oplot, ReadTimes(FileNum*TReads+1:(FileNum+1)*TReads-1), $ PixMDarkHstry(FileNum*TReads+1:(FileNum+1)*TReads-1), $ color=fsc_color('black') oplot, ReadTimes(FileNum*TReads+1:(FileNum+1)*TReads-1), $ Intercepts(FileNum)+ Slopes(FileNum)*$ (Indgen(TReads-1)+1)*KeyStr.ExpTime,$ color=fsc_color('blue') xyouts, ReadTimes(FileNum*TReads), $ Min(PixMDarkHstry(PlotInds))+400,$ Strtrim(Slopes(FileNum)), color=fsc_color('black') EndFor If KeyStr.TvFlag eq 1 and FirstPixel then stop FirstPixel = 0 EndIf EndFor EndFor End ;***************************************************Average Slope per Radii 1: Begin If KeyStr.TvFlag eq 2 then $ device, filename = KeyStr.PlotDir+'M16_'+KeyStr.ThisFilter+$ '_AveragedSlopesRadialPersistence.ps' NumDithers = KeyStr.NumDithers ;Loop through the stars and do the average slope values For StarNum = 0, NumStars-1 do begin print, 'Working on star : '+Strtrim(StarNum,2) XSigs = Fltarr(NumStars) & YSigs = Fltarr(NumStars) XStart = Round(XStar(StarNum)) - 25 > 0 XStop = Round(XStar(StarNum)) + 25 < KeyStr.Naxis1 YStart = Round(YStar(StarNum)) - 25 > 0 YStop = Round(YStar(StarNum)) + 25 < KeyStr.Naxis1 Centroid = Find_Centroid_Sat(SatMask, XStart, XStop, YStart, YStop, $ IsStar, SatRad, TvFlag = KeyStr.TvFlag) XSigs(StarNum) = Centroid(2) & YSigs(StarNum) = Centroid(3) ;Now the centroid is confirmed and it's a star, look at the radii SatRegSize = 101 & RadMax = sqrt(2*(SatRegSize/2)^2) RadMax = 50 ;Check if we'll stay in bounds If Centroid(0) + SatRegSize/2 gt KeyStr.Naxis1 or $ Centroid(0) - SatRegSize/2 lt 0 or $ Centroid(1) + SatRegSize/2 gt KeyStr.Naxis2 or $ Centroid(1) - SatRegSize/2 lt 0 then goto, NewStar RelXCoords = Indgen(SatRegSize,SatRegSize) mod SatRegSize RelXCoords = RelXCoords - Mean(RelXCoords) RelYCoords = reverse(Indgen(SatRegSize,SatRegSize),2) / SatRegSize RelYCoords = RelYCoords - Mean(RelYCoords) RsqCoords = RelXCoords^2+RelYCoords^2 ; Make an array to hold the ratios and loop through each ; radius for this star, averaging the ratio of the slopes to the last AvgSlopeRatios = Fltarr(RadMax,KeyStr.NumDithers) For Rad = 0, RadMax - 1 do begin print, 'At a radius of : ' +Strtrim(Rad,2) RelPixInRad = where(sqrt(RsqCoords) ge Rad $ and sqrt(RsqCoords) lt Rad+1) AbsPixInds = RelXCoords(RelPixInRad) + Centroid(0) + $ (RelYCoords(RelPixInRad) + Centroid(1))* KeyStr.Naxis1 SlopeRatios = [[Fltarr(NumDithers)]] For PixNum = 0, N_Elements(AbsPixInds) - 1 do begin XInd = long(AbsPixInds(PixNum)) mod KeyStr.Naxis1 YInd = long(AbsPixInds(PixNum)) / KeyStr.Naxis2 If Masks.BadMask(XInd, YInd) ne 1 then begin Return_Pixel_History, XInd, YInd, Times, PixHstry, $ PixMDarkHstry, ReadTimes, Slopes, Intercepts SlopeRatios = [[SlopeRatios],$ [Slopes/Slopes(N_Elements(Slopes)-1)]] EndIf EndFor ;Remove the zero element array at the beginning If (Size(SlopeRatios))[0] gt 1 then $ SlopeRatios = SlopeRatios(*,1:(Size(SlopeRatios))[2]-1) SizeSlopeRatios = Size(SlopeRatios) If SizeSlopeRatios(0) eq 1 then begin AvgSlopeRatios(Rad,*) = SlopeRatios EndIf else begin For SlopeNum = 0, N_Elements(Slopes)-1 do begin djs_iterstat, SlopeRatios(SlopeNum,*), mean = m, $ sigma = sig, mask = ma GoodSlopeRats = where(Ma ne 0, NumGoodSlopeRats) If NumGoodSlopeRats gt 0 then begin AvgSlopeRatios(Rad, SlopeNum) = $ Mean(SlopeRatios(SlopeNum, GoodSlopeRats)) EndIf Else Begin AvgSlopeRatios(Rad, SlopeNum) = $ Mean(SlopeRatios(SlopeNum, *)) EndElse EndFor EndElse EndFor If KeyStr.TvFlag gt 0 then begin ;First plot the dither after saturation for good y limits Plot, AvgSlopeRatios(*,1), $ title = 'Ratio of Slopes to Slope in last Dither ' + $ ' (m[r,d] : m[r,8])' ,$ xtitle = 'Radius from center of star', $ ytitle = 'Slope Ratio', $ background = fsc_color('white'), color= fsc_color('black') OffsetTimes = Times-Times(0) Items = ['Dither 0: 0s', $ 'Dither 1: ' + Strtrim(OffsetTimes(1),2)+' s'] Colors = [ColorNum.ColorNums(1), ColorNum.ColorNums(1)] Psym = [0,0] LineStyles = [1,0] TextColors = [fsc_color('black'),fsc_color('black')] oplot, AvgSlopeRatios(*,0), linestyle = 1, $ color = fsc_color('black') For DitherNum = 2, NumDithers - 1 do begin OPlot, AvgSlopeRatios(*,DitherNum), $ Color = ColorNum.ColorNums(DitherNum) Items = [Items, 'Dither '+Strtrim(DitherNum,2) + ': ' + $ Strtrim(OffsetTimes(DitherNum),2) +' s'] Colors = [Colors, ColorNum.ColorNums(DitherNum)] Psym = [Psym, 0] TextColors = [TextColors, fsc_color('black')] LineStyles = [LineStyles, 0] EndFor legend,items,colors=colors,psym=psym,textcolors=textcolors,$ linestyle = linestyles, position = [0.6,0.85],/normal xyouts, 0.6,0.9, 'Star at : ' + Strtrim(XInd,2) + $ ',' + Strtrim(YInd,2), $ /normal, color = fsc_color('black') xyouts, 0.6,0.85, 'Radius of Saturation : ' + Strtrim(SatRad,2), $ /normal, color=fsc_color('black') EndIF EndFor End 2: Begin ;***************************************************Plot Ratios of Slopes ;Pixels that should be in the region of the star are now flagged. ;Use an And with the original mask to get rid of doubly flagged stars SatPixels = Bytarr(KeyStr.Naxis1, KeyStr.Naxis2) For StarNo = 0, N_Elements(XStar) - 1 do begin MinX = floor(XStar(StarNo)) - 20 > 0 MaxX = floor(XStar(StarNo)) + 20 < KeyStr.Naxis1 MinY = floor(YStar(StarNo)) - 20 > 0 MaxY = floor(YStar(StarNo)) + 20 < KeyStr.Naxis2 SatPixels(MinX:MaxX, MinY:MaxY) = 1 EndFor SatStars = SatMask*SatPixels ; Now make arrays for each of these pixels SatStarInds = where(SatStars ne 0, NumPersist) SatStarX = SatStarInds mod KeyStr.Naxis1 SatStarY = SatStarInds / KeyStr.Naxis1 NumSatReads = SatMask(SatStarInds) PixHstry = UIntArr(NumPersist, NumRawFiles*TReads) PixHstryMDark= LonArr(NumPersist, NumRawFiles*TReads) ;Get the value of the saturated pixels for all of the dithers ;Read a pixel from file at a time (seems to be faster than reading ;entire datacube for less than 10,000 pixels or so For FileNum = 0, NumRawFiles - 1 do begin Print, 'Getting Pixel History from File Number ' +Strtrim(FileNum,2) ;Loop through all of the pixels in this file For SatPixs = 0, NumPersist-1 do begin Fits_Read_DataCube, RawFileNames(FileNum), Im, $ Offset = KeyStr.Offset, $ XStart = SatStarX(SatPixs), XStop = SatStarX(SatPixs), $ YStart = SatStarY(SatPixs), YStop = SatStarY(SatPixs), $ ZStart = 0, ZStop = TReads-1 PixHstry(SatPixs, FileNum*TReads:$ (FileNum+1)*TReads-1)=$ Im(0,0,*) EndFor EndFor ;This is the number of reads in a TReads exposure that ;the pixel was saturated SatReadsCount = 4 If KeyStr.TvFlag eq 2 then device, filename = KeyStr.PlotDir+$ KeyStr.ObjectName+'_'+Strtrim(TReads,2)+$ '_'+KeyStr.ThisFilter+'_PersistenceOfPixelsSaturatedIn_'+$ Strtrim(SatReadsCount,2)+'_Reads.ps' For ReadNum = 0, TReads -1 do begin Plot, [0], [0], /nodata, xrange=[0, NumRawFiles], $ yrange=[-500,3000+ReadNum*100] ;Plot the difference between the biases and the last bias BiasDiffs = [LonArr(NumRawFiles-1)] For SatReads = SatReadsCount-1, SatReadsCount -1 do begin ;Examine the case where SatReads of the reads were saturated ThisSatNumInds = where(NumSatReads eq SatReads, ThisSatNum ) For SatPixs = 0, ThisSatNum-1 do begin ;Find the number of pixels where the number of saturated reads matches SatPixInd = ThisSatNumInds(SatPixs) ;Get the Nth Read for all exposures in the dither ThisReadHist = $ PixHstry(SatPixInd, $ Bindgen(NumRawFiles)*TReads+ReadNum) ;Check to make sure bright object is not found in later pixels djs_iterstat, ThisReadHist, sigrej = 4, sigma = readsig, $ mean = readmean, mask = readmask IllPix = where(readmask(1:NumRawFiles-1) eq 0, NumIll) If NumIll eq 0 then begin ThisBiasDiffs = long(ThisReadHist(1:NumRawFiles-1)) - $ long(ThisReadHist(NumRawFiles-1)) BiasDiffs = [[BiasDiffs], [ThisBiasDiffs]] oplot, Bindgen(NumRawFiles-1)+1, ThisBiasDiffs, psym= 2 EndIf EndFor EndFor EndFor End EndCase If KeyStr.TvFlag eq 2 then device, /close stop End