2009.12/1
前回までで81枚分imcombineした偶奇それぞれの画像(mcomb1124odd.fits,mcomb1124even.fits)が出来たが、これは本当に単純に81枚を足し合わせただけ。もっと画像の精度をあげるためには1枚1枚のflameについて、noiseの大きいflameからの寄与は小さくnoiseの小さいflameからの寄与を大きくなるようにimcombineしないといけない。それを実現するためには統計的に、flame毎に1/σ^2をかけたものをimcombineすればいい。(σ=標準偏差)
また次回以降等級を測ったりするが、SExtractorで出してくれるcount=fluxから等級に変換するためにはzero点が必要。zero点決めのためにはいわゆる標準星の情報が必要になるので今日はとりあえず標準星画像のflat-fieldingのみ行う。
☆其の一 各flame毎にσを求め,そのσを使って重み付けして足し合わせ
◎各flameの標準偏差求める
・各フレーム毎のσを求めるcl scriptを作る
cp maskodd.cl sigma1201odd.cl
emacs sigma1201odd.cl &
右のように書き換えて保存→sigma1201odd.cl
偶数チップ分も同様にσ求めるclスクリプト組む
→sigma1201even.cl
・clスクリプトの実行
(imfill使えるようにパッケージstsdasに入っとく)
cl> stsdas
実行
st> cl< sigma1201odd.cl
→sigdisshift38???.fits(位置合わせした80枚の天体flame。maskかかったpixelのvalue=-1000000。別にいらないが)、各flameごとの重み(1/σ^2)のリストsigmaodd.datが出来た！
偶数チップ分も同様に
st> cl< sigma1201even.cl
◎この重み1/σ^2を使って重みづけして位置合わせした81枚をimcombine
st> epar imcombine
I R A F
Image Reduction and Analysis Facility
PACKAGE = immatch
TASK = imcombine
input = @disshift1117odd List of images to combine
output = omomi1201odd.fits List of output images
(headers= ) List of header files (optional)
(bpmasks= ) List of bad pixel masks (optional)
(rejmask= ) List of rejection masks (optional)
(nrejmas= ) List of number rejected masks (optional)
(expmask= expomomi1201odd ) List of exposure masks (optional) #exposure mask(各pixelが何秒or何枚分imcombineされたか＝valueになってる。結果をdisplayさせて見ればわかる)をつくっとく#
(sigmas = sigmamapodd ) List of sigma images (optional) #各pixelごとのカウント値のばらつきを表示。ポアソン分布に従えばシグナル大きいpixelほどノイズも大きいので、天体写ってるpixelほどsigma imageのvalueも大きい。#
(logfile= STDOUT) Log file
(combine= average) Type of combine operation
(reject = sigclip) Type of rejection
(project= no) Project highest dimension of input images?
(outtype= real) Output image pixel datatype
(outlimi= ) Output limits (x1 x2 y1 y2 ...)
(offsets= none) Input image offsets
(masktyp= goodvalue) Mask type #maskがけ忘れずに#
(maskval= 0.) Mask value
(blank = 0.) Value if there are no pixels
(scale = none) Image scaling
(zero = none) Image zero point offset
(weight = @sigmaodd.dat) Image weights #ここに先に作った重みづけの1/σ^2のファイルを設定しとく#
(statsec= ) Image section for computing statistics
(expname= EXPTIME) Image header exposure time keyword #exposure mapの種類を積分時間にしとく#
(lthresh= -1000.) Lower threshold
(hthresh= INDEF) Upper threshold
(nlow = 1) minmax: Number of low pixels to reject
(nhigh = 1) minmax: Number of high pixels to reject
(nkeep = 1) Minimum to keep (pos) or maximum to reject (neg)
(mclip = yes) Use median in sigma clipping algorithms?
(lsigma = 3.) Lower sigma clipping factor
(hsigma = 3.) Upper sigma clipping factor
(rdnoise= 0.) ccdclip: CCD readout noise (electrons)
(gain = 1.) ccdclip: CCD gain (electrons/DN)
(snoise = 0.) ccdclip: Sensitivity noise (fraction)
(sigscal= 0.1) Tolerance for sigma clipping scaling corrections
(pclip = -0.5) pclip: Percentile clipping parameter
(grow = 0.) Radius (pixels) for neighbor rejection
(mode = ql)
(:go)
・これで重み付けをして足し合わせた画像omomi1201odd.fitsとexposure map=expomomi1201odd.plとsigmaマップが出来たので確認してみる
st> !ds9 &
st> reset stdimage=imt4096
st> display omomi1201odd.fits 1 zr- zs- z1=0 z2=200 →こんなの
st> display mcomb1124odd.fits 2 zr- zs- z1=0 z2=200 →こんなの
・・・ぶっちゃけ重みづけしようとしまいと大して変わらない。
exposure map見てみる
st> display expomomi1201odd.pl 3 →こんなの
z1=0. z2=12150.
真中ほど多数のファイルが重なっている(合計の露光時間長い)ことがよくわかる
次にsigmaマップも見てみる
st> display sigmapodd 4 →こんなの
count高い所=天体の場所ほどsigmaが大きいというポアソン分布の確認になってる
偶数チップ分も同様にimcombine
st> epar imcombine
I R A F
Image Reduction and Analysis Facility
PACKAGE = immatch
TASK = imcombine
input = @disshift1117even List of images to combine
output = omomi1201even.fits List of output images
(headers= ) List of header files (optional)
(bpmasks= ) List of bad pixel masks (optional)
(rejmask= ) List of rejection masks (optional)
(nrejmas= ) List of number rejected masks (optional)
(expmask= expomomi1201even ) List of exposure masks (optional) #exposure mask(各pixelが何秒or何枚分imcombineされたか＝valueになってる。結果をdisplayさせて見ればわかる)をつくっとく#
(sigmas = sigmamapeven ) List of sigma images (optional) #各pixelごとのカウント値のばらつきを表示。ポアソン分布に従えばシグナル大きいpixelほどノイズも大きいので、天体写ってるpixelほどsigma imageのvalueも大きい。#
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・（略）・・・・・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
(weight = @sigmaeven.dat) Image weights #ここに先に作った重みづけの1/σ^2のファイルを設定しとく#
(statsec= ) Image section for computing statistics
(expname= EXPTIME) Image header exposure time keyword #exposure mapの種類を積分時間にしとく#
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・（略）・・・・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
(:go)
こちらも確認させてみると
st> display omomi1201even.fits 5 zs- zr- z1=0 z2=200 →こんなの
st> display mcomb1124even.fits 6 zs- zr- z1=0 z2=200 →こんなの
st> display expomomi1201even.pl 7 →こんなの
以上で偶奇チップ毎に、「flatわり,sky引き,位置合わせ,maskがけ,1/σ^2で重みづけして足し合わせ」の済んだ精度のよい画像が作れた。(omomi1201odd.fits,omomi1201even.fits)
以降はこの二枚を足し合わせる方向に進んでいく。
☆其の二 標準星のリストをつくりflat fieldingする
チップ１とチップ２の画像を足すのは天体の位置合わせをしてimcombineすれば良いだけと思ってしまうが、チップ１とチップ２では感度が違うため同じ天体を見てもcountが違う。つまりチップ１とチップ２のzero点(1count/sが何magに対応するか?)が異なるのである。よって二枚の画像についてzero点を合わせてから足し合わせないといけない。
それを本格的に行うのは次週になるが、今日はまずzero点決めにかかせない標準星の画像の一次処理を、とりあえずflat-fieldingだけ行う
・cl0024_data_selectから標準星のデータだけのリストつくる
grep FS cl0024_data_select | cat -n | awk '$1%2==1{print $2,$3,$4,$5,$6,$7}' > standardodd //元データのflameid,object,filter01,exptime,coaddが列挙されてるリストcl0024_data_selectの中から標準星のデータを選んで(FSがある行だけとってくる)、さらにそのうちの奇数番目だけのリストつくる//
grep FS baka2 | cat -n | awk '$1%2==0{print $2,$3,$4,$5,$6,$7}' > standardeven
fitsファイルが羅列された作業用リストもつくる
awk '{print $1".fits"}' standardodd > FSodd
awk '{print $1".fits"}' standardeven > FSeven
・flat fieldingしたあとの結果リストを用意(中身はflFSMCSA38???.fits)
$ awk '{print "flFSMCSA"substr($1,8,20)}' FSodd > flFSodd
$ awk '{print "flFSMCSA"substr($1,8,20)}' FSeven > flFSeven
・以前作ったflat flame=flat1006odd.fits,flat1006even.fitsを使って標準星をflat-fieldingしてやる
st> imarith @FSodd / flat1006odd.fits @flFSodd //以前作ったflat flameで割ってやる→flat補正された標準星ファイル(flFSMCSA38???.fits)が80枚入ったリストflFSoddの出来上がり//
偶数チップ分も同様に
st> imarith @FSeven / flat1006even.fits @flFSeven
結果を一応確認してみる
cl> display flFSMCSA38615.fits 1 zr- zs- z1=0 z2=3000
z1=0. z2=3000.
cl> display MCSA00038615.fits 2 zr- zs- z1=0 z2=3000
z1=0. z2=3000.
よく分からないがflFSMCSA38???.fitsはMCSA00038???.fitsよりもムラがないようにも見える。
