##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009 ##PlotRadialProfiles.py ## ##PURPOSE: ##To plot a series of radial profiles to try and understand the persistence ## ##INPUTS: ## FitsFileName: str ## The full path to the fits file that will be plotted ## ##KEYWORDS: ## UseRef: int ## 0-Don't subtract reference ## 1-Do subtract ## PlotMode: int ## 0 - Plot ramps individually, clearing plot each time ## 1 - Plot all the ramps over each other ## DarkSub: int ## 0 - Don't subtract a dark ## 1 - Subtract a median dark with the same NReads ## DoCentroid: int ## 0 - Don't center image around intenisty peak ## 1 - Center the image around intensity ## XCen: int ## X center coordinate for region of interest ## YCen: int ## Y center coordinate for region of interest ## BoxSize: int ## Size of box centered around (XCen, YCen) ## RadMode: int ## 0 - Take average values at given radius ## 1 - Return values for all radii; no averaging ## Slope: int ## 0 - The datacube is a genuine datacube with NAxis3 up the ramp reads ## 1 - The datacube contains a slopefit; the first slice is the image ## AttFit: int ## 0 - Don't try to fit the radial function ## 1 - Try to fit the radial function with a function of the form ## ## A*Cos(B*r^C)*Exp(-r^2/D) ## Mutli: int ## 0 : do only one plot (in this case, FitsFileName, XCen, and YCen should be single values) ## 1 : do plots from multiple files (in this case, FitsFileName, XCen, and YCen should be arrays) ## ##CALLING SEQUENCE: ## run PlotRadialProfiles.py [FitsFileName] ## ##EXAMPLES: ##________ ##SLOPE IMAGE: ## run PlotRadialProfiles.py '/nfs/slac/g/ki/ki03/lances/H1RG-022/ASIC/Reduced/07Nov19/SAO5417/SAO5417_I_NegPerDead_H1RG_SIPIN_2500_Reads_Nov19_2007_22_15_46_I_RPS_SlopeFlatFielded.fits' --YCen=254 --XCen=911 --Slope=1 --BoxSize=100 ##MULTIPLE FILES: ## run PlotRadialProfiles.py '/nfs/slac/g/ki/ki04/lances/H2RG-32-147/ASIC/07Dec19/NGC956/Saturation_45_Files.lst' --Multi=1 ##//////////////////////////////////////////////////////////////////////// ##SETUP execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py') execfile('/afs/slac/u/ki/lances/python/SPIEPlotSettings.py') from HxRG_Class import * from ReturnCentroid import * from ReturnRadialAvg import * from ReturnRadialProfiles import * from ReadCol import * from Return_Persistence_Quantities import * import TableIO import csv import pdb import time #Get the filename and options from command line FitsFileName = sys.argv[1] ##PARSE KEYWORDS keywords = ['UseRef=','DarkSub=','PlotMode=',\ 'DoCentroid=','XCen=', 'YCen=', 'BoxSize=', \ 'RadMode=','AttFit=', 'Mode=', 'TvFlag=', 'Slope=', \ 'Multi=', 'SlopeMode='] UseRef = 0 DarkSub = 0 PlotMode = 0 DoCentroid = 1 XCen = 100 YCen = 100 BoxSize = 50 RadMode = 1 AttFit = 0 Mode = 0 Slope = 0 TvFlag = 0 Multi = 1 SlopeMode = 0 opts, extraparams = getopt.getopt(sys.argv[2:],'',keywords) for o,p in opts: if o in ['--UseRef']: UseRef = int(p) elif o in ['--DarkSub']: DarkSub = int(p) elif o in ['--PlotMode']: PlotMode= int(p) elif o in ['--DoCentroid']: DoCentroid=int(p) elif o in ['--BoxSize']: BoxSize = int(p) elif o in ['--RadMode']: RadMode = int(p) elif o in ['--XCen']: XCen = int(p) elif o in ['--YCen']: YCen = int(p) elif o in ['--AttFit']: AttFit = int(p) elif o in ['--Mode']: Mode = int(p) elif o in ['--Slope']: Slope = int(p) elif o in ['--TvFlag']: TvFlag = int(p) elif o in ['--Multi']: Multi = int(p) elif o in ['--SlopeMode']: SlopeMode = int(p) #############################################################DO JUST ONE PLOT FROM ONE FILE if Multi == 0: ##GET FILES HxRG = HxRG_C(FitsFileName=FitsFileName) XCen = int(XCen[0]) YCen = int(YCen[0]) if SlopeMode == 1: Radii, IntArr, Read=ReturnRadialProfiles(FitsFileName, Mode, TvFlag=1,\ XCen=XCen, YCen=YCen, Slope=Slope, BoxSize=BoxSize,\ DoCentroid=DoCentroid) else: #Determine coordinates XStart=XCen-BoxSize/2 ; XStop=XCen+BoxSize/2 YStart=YCen-BoxSize/2 ; YStop=YCen+BoxSize/2 PlotDir = HxRG.DataDir+HxRG.DetStr+'/Plots/NegPersistence/' ##FIND THE MATCHING DARK if DarkSub == 1: HxRG.Get_Median_Dark(HxRG.NReads) DarkHDU = pyfits.open(HxRG.DarkName) ##USE THESE SYMBOLS AND COLORS PSyms = ['r+','b+','g+','m+','c+','FF+'] PlotHandles = [] LegendLabs = [] LegendSyms = [] FiltersArr = [] FiltersStr = '' ##DEFINE SOME FITTING FUNCTIONS ExpFunc = 1 if ExpFunc == 1: FitFunc = lambda p, x: p[0]*cos(p[1]*x**p[3])*exp(-x**2/p[2]) else: FitFunc = lambda p, x: p[0]*cos(p[1]*x**p[3])*x**(p[2]) ErrFunc = lambda p, x, y: FitFunc(p,x)-y ##GO THROUGH THE FILES IN SEQUENCE FitsHDU = pyfits.open(FitsFileName) HxRG = HxRG_C(FitsFileName=FitsFileName) HxRG.Get_Raw_Header(FitsFileName) #Find the center of the intensity first for Radial Purposes #Prepare the centroid image axes mplot.figure(0,figsize=[3,3]) mplot.clf() if DoCentroid == 1: Im0 = FitsHDU[0].section[0,:,:] Im1 = FitsHDU[0].section[HxRG.NAxis3-1,:,:] if UseRef == 1: Im0 = ReturnRSFrame(Im0,Im0,HxRG.CFac,BadCols=HxRG.BadCols, HxRG=HxRG) Im1 = ReturnRSFrame(Im1,Im0,HxRG.CFac,BadCols=HxRG.BadCols, HxRG=HxRG) if DarkSub == 1: Dark1 = DarkHDU[0].section[1,:,:] Im1 = Im1-Dark1 Im1 = Im1-Im0 XStart, XStop, YStart, YStop = ReturnCentroid(Im1, XStart, XStop,\ YStart, YStop, TvFlag=1) Read = Im1[YStart:YStop, XStart:XStop] mplot.figure(1,figsize=[3,3]) mplot.clf() Radii, IntArr = ReturnRadialAvg(Read,Mode=1) mplot.plot(Radii,IntArr,'r+') mplot.figure(2) #Form an array to hold all of the mins, maxes, etc. MaxArrs = zeros([HxRG.NAxis3-1],dtype=float32) MinArrs = zeros([HxRG.NAxis3-1],dtype=float32) FZeros = zeros([HxRG.NAxis3-1],dtype=float32) SZeros = zeros([HxRG.NAxis3-1],dtype=float32) MinRads = zeros([HxRG.NAxis3-1],dtype=float32) MaxRads = zeros([HxRG.NAxis3-1],dtype=float32) Sums = zeros([HxRG.NAxis3-1],dtype=float32) Avgs = zeros([HxRG.NAxis3-1],dtype=float32) #Loop through all the reads and fill the array for ReadNum in arange(HxRG.NAxis3): ########################################################### #Get Bias for First Read; Subtract Bias for all others if ReadNum == 0 : #Get the bias frame Im0 = FitsHDU[0].section[0,:,:] if UseRef == 1: Im0 = ReturnRSFrame(Im0,Im0,HxRG.CFac,BadCols=HxRG.BadCols, HxRG=HxRG) if DarkSub == 1: Im0 = Im0 Read0 = Im0[YStart:YStop, XStart:XStop] IntArrs = zeros([Read0.shape[0], Read0.shape[1], HxRG.NAxis3],\ dtype=float32) ############################################################ #Return the radial plot if ReadNum > 0 : Im = FitsHDU[0].section[ReadNum,:,:] if UseRef == 1: Im = ReturnRSFrame(Im,Im0,HxRG.CFac,BadCols=HxRG.BadCols, HxRG=HxRG) if DarkSub == 1: Dark = DarkHDU[0].section[ReadNum,:,:] Im = Im-Dark Read = Im[YStart:YStop, XStart:XStop] Read = Read-Read0 #Add the radial plot and get the true min/max'es if RadMode == 1: Radii, IntArr = ReturnRadialAvg(Read,Mode=RadMode) else: IntArr3 = ReturnRadialAvg(Read,Mode=RadMode) Radii= IntArr3[:,0] IntArr = IntArr3[:,1] if TvFlag == 1: mplot.plot(Radii, IntArr) Rad1d = Radii.copy() Rad1d.shape = Rad1d.size IntArr1d = IntArr.copy() IntArr1d.shape = IntArr1d.size MinInt = IntArr1d.min() MaxInt = IntArr1d.max() MinIntInd = (where(IntArr1d == MinInt))[0] MaxIntInd = (where(IntArr1d == MaxInt))[0] MinRad = Rad1d[MinIntInd[0]] MaxRad = Rad1d[MaxIntInd[0]] SumInt = IntArr.sum() ##Fit the Radii if AttFit == 1: p0=[float(MaxInt), 0.3, 40, .5] p1, success = optimize.leastsq(ErrFunc, p0[:], \ args = (Rad1d, IntArr1d),full_output=0,\ maxfev=10000000.) print 'P1:' +str(p1[0])+' , '+str(p1[1])+' , '+str(p1[2]) + ' , '+\ str(p1[3]) if ExpFunc == 1: Fit = p1[0]*cos(p1[1]*Rad1d**p1[3])*exp(-Rad1d**2/p1[2]) else: Fit = p1[0]*cos(p1[1]*Rad1d)*Rad1d**(p1[2]) MinInt = Fit.min() MaxInt = Fit.max() MinIntInd = (where(Fit == MinInt))[0] MaxIntInd = (where(Fit == MaxInt))[0] MinRad = Rad1d[MinIntInd[0]] MaxRad = Rad1d[MaxIntInd[0]] ##Find the zero crossings IntProd = IntArr1d[0:len(IntArr1d)-1]*IntArr1d[1:len(IntArr1d)] Criterium = (IntProd < 0) | (IntProd == 0) ZeroCrossings = where(Criterium)[0] FirstZRad = Rad1d[ZeroCrossings[0]] SecondZRad = Rad1d[ZeroCrossings[1]] mplot.hold(False) mplot.plot(Rad1d,IntArr1d,'r+') mplot.hold(True) if AttFit == 1: mplot.plot(Rad1d,Fit) time.sleep(2) #Fill the arrays MaxArrs[ReadNum-1] = MaxInt MinArrs[ReadNum-1] = MinInt FZeros[ReadNum-1] = FirstZRad SZeros[ReadNum-1] = SecondZRad MinRads[ReadNum-1] = MinRad MaxRads[ReadNum-1] = MaxRad print 'Read: ' +str(ReadNum)+' , '+'Sum: '+str(SumInt) print 'MinRad: ' +str(MinRad) +' , '+'Min: '+str(MinInt) print 'MaxRad: ' +str(MaxRad) +' , '+'Max: '+str(MaxInt) print 'First Zero: ' +str(FirstZRad) print 'Second Zero:' +str(SecondZRad) else: MaxArrs[ReadNum-1] = -1 MinArrs[ReadNum-1] = -1 FZeros[ReadNum-1] = -1 SZeros[ReadNum-1] = -1 MinRads[ReadNum-1] = -1 MaxRads[ReadNum-1] = -1 if RadMode == 1: IntArrs[:,:,ReadNum]=IntArr else: IntArrs[:,:,ReadNum]=Read Sums[ReadNum-1] = Read.sum() Avgs[ReadNum-1] = Read.mean() ###################################################################DO TWO PLOTS FROM DIFFERENT FILES else: FitsFileName, XCen, YCen = ReadDitherList(FitsFileName) HxRG1 = HxRG_C(FitsFileName=FitsFileName[0]) HxRG2 = HxRG_C(FitsFileName=FitsFileName[1]) XCen1 = int(XCen[0]) ; XCen2 = int(XCen[1]) YCen1 = int(YCen[0]) ; YCen2 = int(YCen[1]) if SlopeMode == 1: Radii, IntArr, Read=ReturnRadialProfiles(FitsFileName, Mode, TvFlag=1,\ XCen=XCen, YCen=YCen, Slope=Slope, BoxSize=BoxSize,\ DoCentroid=DoCentroid) else: #Determine coordinates XStart1=XCen1-BoxSize/2 ; XStop1=XCen1+BoxSize/2 YStart1=YCen1-BoxSize/2 ; YStop1=YCen1+BoxSize/2 XStart2=XCen2-BoxSize/2 ; XStop2=XCen2+BoxSize/2 YStart2=YCen2-BoxSize/2 ; YStop2=YCen2+BoxSize/2 PlotDir = HxRG1.DataDir+HxRG1.DetStr+'/Plots/NegPersistence/' ##FIND THE MATCHING DARK if DarkSub == 1: HxRG1.Get_Median_Dark(HxRG1.NReads) DarkHDU1 = pyfits.open(HxRG1.DarkName) ##USE THESE SYMBOLS AND COLORS PSyms = ['r+','b+','g+','m+','c+','FF+'] PlotHandles = [] LegendLabs = [] LegendSyms = [] FiltersArr = [] FiltersStr = '' ##DEFINE SOME FITTING FUNCTIONS ExpFunc = 1 if ExpFunc == 1: FitFunc = lambda p, x: p[0]*cos(p[1]*x**p[3])*exp(-x**2/p[2]) else: FitFunc = lambda p, x: p[0]*cos(p[1]*x**p[3])*x**(p[2]) ErrFunc = lambda p, x, y: FitFunc(p,x)-y ##GO THROUGH THE FILE1 FitsHDU1 = pyfits.open(FitsFileName[0]) HxRG1 = HxRG_C(FitsFileName=FitsFileName[0]) HxRG1.Get_Raw_Header(FitsFileName[0]) ##GO THROUGH THE FILE1 FitsHDU2 = pyfits.open(FitsFileName[1]) HxRG2 = HxRG_C(FitsFileName=FitsFileName[1]) HxRG2.Get_Raw_Header(FitsFileName[1]) #Find the center of the intensity first for Radial Purposes #Prepare the centroid image axes mplot.figure(0,figsize=[3,3]) mplot.clf() if DoCentroid == 1: Im0_1 = FitsHDU1[0].section[0,:,:] Im1_1 = FitsHDU1[0].section[HxRG1.NAxis3-1,:,:] if UseRef == 1: Im0_1 = ReturnRSFrame(Im0_1,Im0_1,HxRG1.CFac,BadCols=HxRG1.BadCols, HxRG1=HxRG) Im1_1 = ReturnRSFrame(Im1_1,Im0_1,HxRG1.CFac,BadCols=HxRG1.BadCols, HxRG1=HxRG) if DarkSub == 1: Dark1_1 = DarkHDU1[0].section[1,:,:] Im1_1 = Im1-Dark1 Im1_1 = Im1_1-Im0_1 XStart1, XStop1, YStart1, YStop1 = ReturnCentroid(Im1_1, XStart1, XStop1,\ YStart1, YStop1, TvFlag=1) Read_1 = Im1_1[YStart1:YStop1, XStart1:XStop1] mplot.figure(1,figsize=[3,3]) mplot.clf() Radii1, IntArr1 = ReturnRadialAvg(Read_1,Mode=1) mplot.plot(Radii1,IntArr1,'r+') mplot.figure(2) #Loop through all the reads and fill the array for ReadNum in arange(HxRG1.NAxis3): ########################################################### #Get Bias for First Read; Subtract Bias for all others if ReadNum == 0 : #Get the bias frame Im0_1 = FitsHDU1[0].section[0,:,:] Im0_2 = FitsHDU2[0].section[0,:,:] if UseRef == 1: Im0_1 = ReturnRSFrame(Im0_1,Im0_1,HxRG1.CFac,BadCols=HxRG1.BadCols, HxRG=HxRG1) Im0_2 = ReturnRSFrame(Im0_2,Im0_2,HxRG1.CFac,BadCols=HxRG1.BadCols, HxRG=HxRG1) if DarkSub == 1: Im0_1 = Im0_1 Im0_2 = Im0_2 Read0_1 = Im0_1[YStart1:YStop1, XStart1:XStop1] Read0_2 = Im0_2[YStart2:YStop2, XStart2:XStop2] IntArrs1 = zeros([Read0_1.shape[0], Read0_1.shape[1], HxRG1.NAxis3], dtype=float32) IntArrs2 = zeros([Read0_2.shape[0], Read0_2.shape[1], HxRG2.NAxis3], dtype=float32) ############################################################ #Return the radial plot if ReadNum > 0 : Im_1 = FitsHDU1[0].section[ReadNum,:,:] Im_2 = FitsHDU2[0].section[ReadNum,:,:] if UseRef == 1: Im_1 = ReturnRSFrame(Im1,Im0_1,HxRG1.CFac,BadCols=HxRG1.BadCols, HxRG=HxRG1) Im_2 = ReturnRSFrame(Im2,Im0_2,HxRG2.CFac,BadCols=HxRG2.BadCols, HxRG=HxRG2) if DarkSub == 1: Dark1 = Dark1HDU[0].section[ReadNum,:,:] Im_1 = Im_1-Dark1 Dark2 = Dark2HDU[0].section[ReadNum,:,:] Im_2 = Im_2-Dark2 Read1 = Im_1[YStart1:YStop1, XStart1:XStop1] Read1 = Read1-Read0_1 Read2 = Im_2[YStart2:YStop2, XStart2:XStop2] Read2 = Read2-Read0_2 #Add the radial plot and get the true min/max'es if RadMode == 1: Radii1, IntArr1 = ReturnRadialAvg(Read1,Mode=RadMode) Radii2, IntArr2 = ReturnRadialAvg(Read2,Mode=RadMode) else: IntArr3_1 = ReturnRadialAvg(Read1,Mode=RadMode) IntArr3_2 = ReturnRadialAvg(Read2,Mode=RadMode) Radii1 = IntArr3_1[:,0] IntArr1 = IntArr3_1[:,1] Radii2 = IntArr3_2[:,0] IntArr2 = IntArr3_2[:,1] if TvFlag == 1: mplot.clf() mplot.hold(True) mplot.plot(Radii1, IntArr1, 'ro') mplot.plot(Radii2, IntArr2, 'bo') pdb.set_trace()