1次処理メモ (可視高分散分光データ)


UNDER CONSTRUCTION・・(誤りがあるかもしれない)＊すばる望遠鏡HDSや岡山188cm望遠鏡HIDESなどの分光器のデータ解析が念頭。＊光赤天連のwikiや 森谷さんのページを参考にさせて頂いた、自分用のメモ＊詳しいパラメータ等は、Subaru Data Reduction CookbooksやSubaru/HDS webpageを参照。＊Flux補正等は未勉強・・

＊ irafコマンドメモ(splotなど)も参考。

観測で取得するフレーム

Objectフレーム観測したい天体を撮った生データBIASフレーム読み出しノイズを評価するために、露出時間0秒で読み出しのみを行ったデータ。Darkフレーム長時間露出をかけると、光が当たっていないのに勝手に励起してくる熱電子がカウントとして沢山載ってくる場合がある。この影響を除くために、シャッターを閉じた状態でObjectフレームと同じ露出時間を使って得た画像。Subaru/HDSやOkayama/HIDESでは、このDark成分のカウントは十分小さいことが分かっているので、撮らない。Comparisonフレーム波長較正に用いる。波長の基準。Th-Arランプなどを用いる。(Subaru HDSは、Th-Arの輝線同定)FLATフレーム各pixel毎の感度ムラの補正用に、検出器に対して一様光源を当てて得られた画像。Image Slicerを観測で用いている場合、slitモードでも取得しておくとよいようだ。


解析手順

1.Overscan

Subaru/HDSのケース・irafの「noao/imred/echelle/」へ入る。・overscan用のスクリプト(overscan.cl)を観測所HPからダウンロードして来て、処理をしているディレクトリに置く。・「ls *.fits > over-in.list」とかで入力の際のリストファイルを作成すると非常に便利。(以下もリストファイルを作って作業する)・task overscan=overscan.cl・overscan *.fits[0] os*.fits
OAO/HIDESのケース・irafの「noao/imred/echelle/bias」へ入る。・overscan用の簡易スクリプト(quick_overscan.cl：神戸さん提供)を処理をしているディレクトリに置く。・task quick_overscan=quick_overscan.cl・quick_overscan h**** ovh**** (****はd0267777など) (これで一度にbiasまで引ける)・irafの「noao/imred/echelle/」へ戻る。

この後、整理の為、blueとredなどで、ディレクトリを分けて解析を進めると良い。ちなみに、「implot *.fits」で、スペクトルの2次元断面見れる(flatカウント確認等)表示範囲の変更は、「:x 1000 2000」(横軸ピクセル)や「:y 0 5000」(縦軸カウント)のようにする。

2.Bias引き(+Dark引き)

・HIDES-Fでは、quick_overscan.clをつかっていれば、overscan処理と同時にbiasを引いてしまっている。・dark成分は、HDSやHIDESではほとんど効かないので、無視してよい。
処理内容・imcombineでbiasフレームを足し合わせ(注：medianをとる)⇒「bias.fits」・imarith ov*.fits - bias.fits bi-ov*.fits ・前半夜後半夜で2回撮っている場合、取りあえず足して処理している。

3.スペクトル抽出(参照の為)

apallを使う星、flatそれぞれ行う(注：星とflatでは幅や光の当たり方が微妙に違うし、幅も違う)。HDSのImage SlicerやHIDES-F(Fiber)を用いている時は、切り出しの幅やピークからどの高さまで取ってくるかに特に注意する。(赤い方のフレームと青い方のフレームで切り出しのパラメータは違う。あまり広げ過ぎると、今度は変なノイズを拾う事も.....)
※HDSでimage slicer使用の場合参照用としては、Slit長さを狭くしたflatを使うcaseも。※「extras」,「bkg」は「no」に、「avglimi」は「yes」としておく。※traceの幅を決める位置は、[line]の数字を調整する※apeture設定時：「d」でaperture削除、「n」でapperture追加。 番号付け直し：1番にしたい所にカーソルを持っていき、「o」のあと、質問に対し「1」と入力

4.flat fielding

・imcombineで足し合わせてフラットフレームを作る (flat.fits)
・apnormalizeでフラットを規格化 (flatn.fits) 3で抽出したスペクトルを、maskingのreferenceにする apresizeで切り出し幅設定：flatの場合の切り出し幅は星より少し狭くしておく。
・規格化したフラット(flatn.fits)で、観測データを割る imarith **.fits / flatn.fits flfi-*.fits

5.背景光除去

apscatter 3.で抽出したスペクトルを参照で使う。 「find,trace,fittrac」はno apresizeでmask幅をちゃんと設定。(apallの時よりは広めにする)
slitモードの場合は、まあまあsky成分を引けるが、Image Slicerやfiberでは難。
Image Slicerやfiberでは背景光のカウントが低い(平均カウント＜ばらつきの大きさ)ので、fittingのパラメータをうまく調整(マイナスを引かない！！) (例)：low_rej:5, high_rej:2
当然だが、天体カウントが高いと、背景光は大きい事が分かる(Image Slicerでも)。

6.宇宙線除去

すばるHDSのページにあるマニュアル(リンク)の31‐33ページを参照
○行う事の要点(もとのフレームをAとする) ・medianフィルタ(注も参照)：AからmedianA作成。medianAは元のAより少しなまっていて、宇宙線はほぼ載っていないimageとなる。 ・a=A/medianA：全体が1カウント程度で、宇宙線ノイズがのったフレームが出来る。 ・aについてlineclean：各行を緩やかにfitし、大きく外れたカウントを消す。 ・imtransでxy軸入れ替えて同作業をし、再びxyを戻す⇒a'作成 ・A'=a'x medianAでAから宇宙線ノイズを除いたA'が出来る。
カウントの引き過ぎに注意する！(注)マニュアルに載っている3×3のmedianより、1×9ののmedianの方がより良く消えた(ただし、処理の最初に500カウント足しておいて、最後に戻すという操作をした。また、1×9なので、linecleanは1方向のみでよい)
自分用メモ上記の手順の実行過程(以下の各ファイルは、本田さん作成の物を若干作り変えたもので、自分の手元にあり。)・背景光除去終了段階のfile listをap-out.listとする。・cr-list.shで、必要なlistファイルを作成・「cl < rejcos-rej1.cl」で実行。「rej1-***」というファイル出来る。

7.1次元化

apall(3で抽出したスペクトルをreferenceにしてトレース)「find,trace,fittrac」はno
HDSのImage SlicerやHIDES-F(Fiber)を用いている時は、切り出しの幅やピークからどの高さまで取ってくるかに特に注意する。 具体的には、「llimit」や「ulimit」,「ylevel」をいじる。
 ＊その時々のSeeing等によって、切り出し幅を変える必要がありそう.... ＊幅を広げ過ぎると、変なガタガタが乗って、カウントが増えても逆にS/Nが悪くなることも。 ＊赤い方のフレームと青い方のフレームでそれぞれ幅が微妙に違う。 ＊暗めの天体と比較星として撮っているA型星などの明るい星のそれぞれで、切り出しのパラメータを幾つかテストしてみる。 ＊「ulimit」や「llimit」、「ylevel」を変える。

8.(オプション) blaze functionの影響を消す(規格化の改善；flatを使用する)

(4で作った)「flat.fits」を「flatn.fits」で割る。(→flfi-flat.fits) それを1次元化。(大まかなblaze functionのデータが得られる。1d-flfi-flat.fits) ＊切りだす際に幅に注意。場合によっては、星とは別のtrace参照データを使う。幅が広すぎると、ガタガタのデータになったりする。 (llmit:-5,ulimit:5,ylevel:0.999とかでよい。特に、Image Slicerでは注意) ⇒⇒slit modeで取ったflatをトレースの参照データとして使うのが、良いかもしれない。 boxcarを使って、各オーダーの平均カウントのデータを作る。(xwindow=分散方向のpix数,ywindow=1, boundar=wrap) (→box-1d-flfi-flat.fits) imarith 1d-flfi-flat.fits / box-1d-flfi-flat.fits nb-1d-flfi-flat.fits
「7.」で1次元化したスペクトル「1d**.fits」を「nb-1d-flfi-flat.fits」で割る。 imarith 1d**.fits / nb-1d-flfi-flat.fits df-1d**.fits ⇒⇒ これによって、そこそこblaze functionが補正され(完全ではないので、重要なラインでどこまでcontinuumが決めやすくなるかは要確認)、 かつカウント数の情報は保持された状態のスペクトルが出来る(S/Nの見積もりもそのままできるので後で便利)。

9.波長較正

まず、Comparisonのデータを一次元化。(referenceは天体。「find,trace,fittrac」に加えて、「recente,resize」もnoとする。)「ecreident」で、Comparisonデータのいずれか一つについて、慎重に波長同定(Th-Ar)
＜波長同定の際に用いるTh-Arのアトラス＞ ・Subaru/HDS link ・OAO/HIDS link 上記のlinkにあるHIDESのAtlasは、約8,500Åまでしかない。それより長波長の分は観測室にあるので、 それをスキャンさせて貰って使うか、HDSのAtlasもうまく活用する。
「ecreident」で、他のComparisonデータについても同定 「ecreident comp2 refe=comp1」
「refspectra」：天体に適用する波長スケールを決定。(sortはutかmjd。groupは空白) 「refspectra objectlist refe=@comparison.list」
「discpor」：波長較正一応comparisonについても行う。

10.複数枚足し合わせ

複数枚足し合わせて、S/Nを上げる。scombine(注group=apertures,combine=sum,reject=sigclip,mclip=yes,lsigma=3,hsigma=1)#また、scombineで、thresholdを設定しておくとよい。(例えば、負のカウントのデータは落す等。)例：「scombine @combine.list combined.fits」

11.Continuumの決定と規格化

一本化までファイル名入力ミスを防ぐために、下記の一連のコマンドを並べたcl scriptを作っておくと便利であった。scombineの設定は、group=images, combine=sum, reject=none,lthresh=0,hthresh=1000000
オーダーの端のS/Nが悪くなるのを防ぐために、以下の方法をとるcontinuum sc*.fits co-sc*.fits (continuumの際の「high_rej」というパラメータは、S/Nの悪いデータでは、3ではなく、2にした方が良いかもしれない。)
scombine sc*.fits cb-sc*.fitsimarith sc*.fits / co-sc*.fits ct-sc*.fits scombine ct-sc*.fits cbct-sc*.fits imarith cb-sc*.fits / cbct-sc*.fits ip*.fits

12.blueとredに分かれていたスペクトルを一つのファイルに。

object-b.fitsとobject-r.fitsの足し合わせのケースで説明。(HIDESのように3枚のCCDの場合も同様)
object-b.fitsとobject-r.fitsの名前が縦に並んだobject-br.listを作る。「scombine @object-br.list object-br.fits」scombineの設定は、group=apertures, combine=averageとする。合わせて、hthreshを設定して、高いカウントはcutした方がよい。

13.太陽中心への補正

fitsのヘッダから観測日、UT、RA、DECを読み取る。「hselect @infits.list DATE-OBS,UT,RA,DEC yes > rvcor.obs」
なお、「rvcor.obs」の中身を2013 01 01 00:00:00.0 00:00:00.0 +00:00:00というような形式にするように注意する。 ※そのままfitsのヘッダを読んでくる、 HDSの場合、日付が2013-01-01のようにハイフンを含んだ形になり、 そのままrvcorをすると変な値になる。「sed -e "s/2013-01-01/2013 01 01/g" rvcor.obs > rvcor-dayrev.obs」
irafのnoao-astutilに移動して、補正速度を求める「rvcor f="rvcor-dayrev.obs"」
なお、rvcorの際に、観測所の設定をする必要がある。 keckの値がirafに入っているので、すばるならkeckにしておけばよい。irafに入っていない観測所の場合以下のようにする。 ・rvcorの設定で、"observa=)_.observatory"とする。 ・observatory taskの設定に観測所の緯度、経度、高度等の情報を入力。
出てくるVHELIOの値を使う。(VHELIOの値だけなら、rvcorのパラメータはirafのdefaultの1900年分点で良さそう)
dopcor @infile @outfile -VHELIO で太陽中心に補正(VHELIOの符号が反転する事に注意！)
dopcorの設定パラメータ epar dopcor input = List of input spectra output = List of output spectra redshift= Redshift or velocity (Km/s) (isveloc= yes) Is the redshift parameter a velocity? (add = no or yesは, スペクトルが1次元データか, 多次元データかで適宜判断) Add to previous dispersion correction? (dispers= yes) Apply dispersion correction? (flux = no) Apply flux correction? (factor = 3.) Flux correction factor (power of 1+z) (apertur= ) List of apertures to correct (verbose= no) Print corrections performed? (mode = ql)


参考：視線速度(別途測定)を補正する場合。「dopcor infile outfile (+の速度)」で青い方向へ動き、「dopcor infile outfile (-の速度)」で赤い方向へ動くので、別に視線速度(RV:redshiftが正、blueshiftが負)を測定していて、実験室系の波長へと持って行く場合、「dopcor infile outfile +RV」とする。
**dopcorのパラメータ「add」は要注意: (以下, help dopcorの記述より)If the add parameter is used and the image uses a "multispec" format where the previous doppler factor is stored separately then the new doppler factor is: znew = (1 + z) * (1 + zold) - 1 = z + zold + z * zold where z is the specified doppler factor, zold is the previous one, and znew is the final doppler factor. If the add parameter is no then the previous correction is replaced by the new correction. Note that for images using a linear or equispec coordinate system the corrections are always additive since a record is not kept of the previous correction. Also any flux correction is made based on the specified doppler correction rather than znew.

参考：HJD(Heliocentric Julian Date)の値をfitsファイルのheaderの情報をもとに算出。irafのnoao-imred-echelleで(あるいはnoao-astutilに移動して)、「setjd」を実行。 設定の例(EXPTIMEを設定した場合、時刻は観測時刻の最初だと判定され、jdは観測中心への補正される)

14.fitsとtxtの変換

noao.onedspecで、fits⇒txt 「wspectext hoge.fits hoge.txt」txt⇒fits 「rspectext hoge.txt hoge.fits」