Skip to content

Commit

Permalink
Merge pull request #2107 from desihub/gains-updates
Browse files Browse the repository at this point in the history 
Minor changes to desi_compute_gains
  • Loading branch information
@sbailey
sbailey authored Aug 30, 2023
2 parents b98b567 + feae0af commit 7dcaa77
Showing 1 changed file with 50 additions and 41 deletions.
91 changes: 50 additions & 41 deletions bin/desi_compute_gain
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,10 +13,10 @@ import matplotlib.pyplot as plt
from desiutil.log import get_logger

def clipped_var(mx,x2) :
nsig=5.

nsig=5.
ovar=0.001

# start with NMAD , less sensitive to outliers
var = ( 1.4826*np.median( np.abs(np.sqrt(x2))) )**2
#print("ini",var,0)
Expand All @@ -28,7 +28,7 @@ def clipped_var(mx,x2) :
break
ovar=var
return np.mean(mx[ok]),var,ok.size,loop


def _parse_sec_keyword(value):
'''
Expand Down Expand Up @@ -56,7 +56,7 @@ parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFo
description="Compute the CCD gains (in e/ADU) using a series of similar images. Images are paired according to their exposure time")
parser.add_argument('-i','--image', type = str, default = None, required = True, nargs = "*",
help = 'path of preprocessed image fits files')
parser.add_argument('--exptime-keyword', type = str, default = "EXPREQ",
parser.add_argument('--exptime-keyword', type = str, default = "EXPTIME",
help = 'change exposure time keyword for pairing images')
parser.add_argument('--max-mean-flux',type = float, default = 8000)
parser.add_argument('--min-pixel-flux',type = float, default = -1000)
Expand All @@ -68,8 +68,8 @@ parser.add_argument('--npy',type = int, default = 1,
help = 'number of adjacent pixels in a column) to add before computing variance')
parser.add_argument('--margin',type = int, default = 200,
help = 'remove margins around CCD')
parser.add_argument('--amplifiers',type = str, default="ABCD", help="amplifiers being studied")
parser.add_argument('--deg',type = int, default=3, help="max degree of polynomial fit")
parser.add_argument('--amplifiers',type = str, default="ABCD", help="amplifiers being studied")
parser.add_argument('--deg',type = int, default=3, help="max degree of polynomial fit")
parser.add_argument('--fix-rdnoise',action='store_true', help="do not refit readnoise")
parser.add_argument('--plot',action="store_true",help="show the fit")
parser.add_argument('--outfile',type = str, default=None, help="save PTC values in ASCII file")
Expand All @@ -84,7 +84,7 @@ camera = None

# loop on preprocessed images
filenames=args.image
for filename in filenames :
for filename in filenames :
hdulist=pyfits.open(filename)
this_camera=hdulist[0].header["CAMERA"].strip()
if camera is None :
Expand All @@ -93,16 +93,25 @@ for filename in filenames :
if this_camera != camera :
log.error("Not the same camera for all images, I find {} and {}".format(camera,this_camera))
sys.exit(1)

exposure_times.append(hdulist[0].header[args.exptime_keyword])

unique_exposure_times = np.unique(exposure_times)
#unique_exposure_times = np.unique(exposure_times)

#log.info("Exposure times = {}".format(unique_exposure_times))


unique_exposure_times = np.array([exposure_times[0]])
threshold=0.2 # sec
for exptime in exposure_times[1:] :
diff=np.min(np.abs(unique_exposure_times-exptime))
if diff>threshold : unique_exposure_times = np.append(unique_exposure_times,exptime)

log.info("Exposure times = {}".format(unique_exposure_times))

pairs = []
for exptime in unique_exposure_times :
ii=np.where(exposure_times==exptime)[0]
ii=np.where(np.abs(exposure_times-exptime)<threshold)[0]
npairs=ii.size//2
for p in range(npairs) :
pairs.append( [ filenames[ii[p*2]],filenames[ii[p*2+1]] ] )
Expand All @@ -123,30 +132,30 @@ gain_tbl['GAIN'] = np.zeros(len(gain_tbl))
gain_tbl['ERRGAIN'] = np.zeros(len(gain_tbl))

for a,amp in enumerate(args.amplifiers) :

ax=[] # mean fluxes
ay=[] # variances for mean fluxes
ardn=[] # variance of readnoise
an=[] # number of data points in bin

nbins=int(args.max_mean_flux/args.bin_size)

for p,pair in enumerate(pairs) :
log.info("pair #{} {} {}".format(p,pair[0],pair[1]))
h1=pyfits.open(pair[0])
h2=pyfits.open(pair[1])
k="CCDSEC%s"%amp

yy,xx=_parse_sec_keyword(h1[0].header[k])

img1=h1[0].data[yy,xx]
img2=h2[0].data[yy,xx]

rdnoise1=h1[0].header["OBSRDN%s"%amp]
rdnoise2=h1[0].header["OBSRDN%s"%amp]
rdnoise=np.sqrt((rdnoise1**2+rdnoise2**2)/2.)
#log.info("rdnoise : {},{} -> {}".format(rdnoise1,rdnoise2,rdnoise))


if args.margin>0 :
margin=args.margin
Expand All @@ -166,21 +175,21 @@ for a,amp in enumerate(args.amplifiers) :

img1[img1>args.max_pixel_flux]=1e40 # kill entry
img2[img2>args.max_pixel_flux]=1e40 # kill entry

if args.npx>1 : # rebinning lines
r=args.npx
n0=img1.shape[0]
n0=img1.shape[0]
n1=(img1.shape[1]//r)*r
img1=img1[:,:n1].reshape(n0,n1//r,r).sum(axis=2)
img2=img2[:,:n1].reshape(n0,n1//r,r).sum(axis=2)
if args.npy>1 : # rebinning columns
r=args.npy
n0=(img1.shape[0]//r)*r
n0=(img1.shape[0]//r)*r
n1=img1.shape[1]
img1=img1[:n0,:].reshape(n0//r,r,n1).sum(axis=1)
img2=img2[:n0,:].reshape(n0//r,r,n1).sum(axis=1)
npix=args.npx*args.npy

if 0 : # simple scalar scaling
ok=(img1>200)&(img1<args.max_pixel_flux)&(img2>200)&(img2<args.max_pixel_flux)
scale=np.exp(np.median(np.log(img1[ok]/img2[ok]))) # better
Expand All @@ -192,7 +201,7 @@ for a,amp in enumerate(args.amplifiers) :
ok=(img1[j]>args.min_pixel_flux*npix)&(img1[j]<args.max_pixel_flux*npix)&(img2[j]>args.min_pixel_flux*npix)&(img2[j]<args.max_pixel_flux*npix)
sum1=np.sum(img1[j][ok])
sum2=np.sum(img2[j][ok])
if sum2>100 :
if sum2>100 :
scale[j]=sum1/sum2 # better
for j in range(img1.shape[0]) :
if np.abs(scale[j]-1)>0.2 :
Expand All @@ -204,48 +213,48 @@ for a,amp in enumerate(args.amplifiers) :
plt.figure("ratio")
j=np.arange(img1.shape[0])
plt.plot(j[scale>0],scale[scale>0],label="pair %d"%p)


ok=(img1>args.min_pixel_flux*npix)&(img1<args.max_pixel_flux*npix)&(img2>args.min_pixel_flux*npix)&(img2<args.max_pixel_flux*npix)
xx=(img1[ok]+img2[ok])/2.
yy=(img1[ok]-img2[ok])**2/2. # we use the difference of images so the var of a single flux = var(diff)/2.
yy=(img1[ok]-img2[ok])**2/2. # we use the difference of images so the var of a single flux = var(diff)/2.
bins=(xx/args.bin_size).astype(int)
ok=(bins>=0)&(bins<nbins)
ubins=np.unique(bins[ok])

x=[]
y=[]
n=[]

for b in ubins :
ok=np.where((xx>=args.bin_size*b)&(xx<args.bin_size*(b+1)))[0]
if ok.size<400 : continue
mean,var,ndata,nloop = clipped_var(xx[ok],yy[ok])
log.debug("flux=%f var=%f n=%d nloop=%d"%(mean,var,ndata,nloop))
log.debug("flux=%f var=%f n=%d nloop=%d"%(mean,var,ndata,nloop))
x.append(mean)
y.append(var)
n.append(ndata)
ax.append(mean)
ay.append(var)
an.append(ndata)
ardn.append(npix*rdnoise**2)


if ofile is not None :
for i in range(len(x)) :
ofile.write("%f %f %d %d %d\n"%(x[i],y[i],n[i],a,p))


# fit the data for this amplifier

ok=(np.array(ay)>0.1) # var
x=np.array(ax)[ok]
y=np.array(ay)[ok]
y0=np.array(ardn)[ok]
n=np.array(an)[ok]
n=np.array(an)[ok]
err=np.sqrt(2./n)*y
xs=np.mean(x)

w=1./(err**2+1.**2)

def myfit(w,x,y,deg) :
Expand Down Expand Up @@ -275,7 +284,7 @@ for a,amp in enumerate(args.amplifiers) :
coef.append(myfit(w,x,y,deg))
gain.append(1./coef[deg-1][1]*xs)
log.info("%s %s deg=%d gain = %4.3f e/ADU"%(camera,amp,deg,gain[deg-1]))

mgain=np.mean(gain)
errgain=np.max(gain)-np.min(gain)
log.info("%s %s gain = $%4.3f \\pm %4.3f$ e/ADU"%(camera,amp,mgain,errgain))
Expand All @@ -287,17 +296,17 @@ for a,amp in enumerate(args.amplifiers) :
#chi2=np.sum((y-y0-mypol(x,coef))**2/err**2)
#ndf=x.size-coef.size
#gain=1./coef[1]*xs

if args.plot :
ms=5
fig=plt.figure("{}-{}".format(camera,amp))
plt.subplot(2,1,1)
plt.errorbar(x,y,err,fmt="o",ms=ms,color="gray")
tx=np.arange(np.min(x),np.max(x),100)
tx=np.arange(np.min(x),np.max(x),100)
my0=np.mean(y0)
for d in range(1,args.deg+1) :
plt.plot(tx,my0+mypol(tx,coef[d-1]),label="deg %d fit, gain=%4.3f"%(d,gain[d-1]))

plt.ylabel("variance(flux)")
plt.legend(loc="upper left",title="%s-%s"%(camera,amp))
plt.subplot(2,1,2)
Expand Down Expand Up @@ -337,6 +346,6 @@ for j in range(len(gain_tbl["AMP"])) :

if args.plot :
plt.show()


#plt.show()

0 comments on commit 7dcaa77

Please sign in to comment.