XMM-Newton衛星のデータ解析の解説

概要

XMM-Newton衛星は欧州の打ち上げたX線天文衛星で、この衛星のデータ解析について易しい日本語で解説した資料が見当たらないので、批判を恐れず自分でまとめてみる事にする。(訂正などあれば優しくご指摘ください。)

下記の設定方法などは2018年7月時点での話になります。bash ユーザーを想定してますので、tcshユーザーは適宜変換ください。コメントにあるように、一部、unscheduleのデータが混じる場合にはコケます(2019年10月発覚, まだ未対応です)。

XMM-Newton 衛星の基礎情報

XMM-Newtonについての基礎情報を取得しよう。全部に目を通すのは不可能だが、斜め読みでもして、どのような情報がどこに散らばってるか把握しておこう。

XMM-Newton のホームページ (ESA)
- ESAのページ。
  - https://www.cosmos.esa.int/web/xmm-newton
  - 情報源泉はこれを読むしかない。
- sas-threads : 役立つツールやコマンド集
  - https://www.cosmos.esa.int/web/xmm-newton/sas-threads
XMM-Newton のホームページ (NASA/GSFC)
- https://heasarc.gsfc.nasa.gov/docs/xmm/xmm_gal_science.html
- NASA/GSFCのページ。ESAのページとの違いはよく分からない。
XMM-Newton のユーザーガイド (ESA)
- https://xmm-tools.cosmos.esa.int/external/xmm_user_support/documentation/sas_usg/USG/node1.html
- ESAのユーザーガイド。かなり細かい。
XMM-Newton の ABC ガイド (NASA/GSFC)
- https://heasarc.gsfc.nasa.gov/docs/xmm/abc/
- NASA/GSFCのユーザーガイド。ESAと比べるとコンパクトで良い。

XMM-Newton 衛星の観測装置

X線望遠鏡が３台あって、２つのミラーは焦点面を回折格子(RGS)とCCD(EPIC)でシェアし、一つはCCD(PN)が専有する。有効面積は、Chandra衛星の約4倍で、角度分解能は Chandra よりも一桁悪くて約15秒角(HPD)である。点源までは小さくない程度に広がっている天体や構造の観測ではバランスのよいX線衛星であろう。通常の解析においては、ミラーが３台で、X線CCDが、EPIC-MOS1, MOS2, PN の３台あることは把握しておこう。

観測装置
- X線CCD : EPIC MOS
  - 一枚のWalter-I型ミラーで集光する光子の半分は回折光子(RGS)に行くので、一台あたりの有効面積は EPIC PN の約半分
  - MOS1, MOS2 の２つのモジュールがある。
- X線CCD : EPIC PN
  - 一枚のWalter-I型ミラーで集光する光子を、全部の占有できるため、一台あたりの有効面積はEPIC MOS の約２倍
- X線回折格子 : RGS
  - 低エネルギー側で高分散分光
- 可視光/紫外線モニター : OM
XMM-Newton の概要についての参考資料
- PLAIN Center News : http://www.isas.jaxa.jp/docs/PLAINnews/115_contents/115_1.html
- - 「XMM-Newton 衛星 X 線スペクトルを用いたペルセウス座銀河団中心部の高温ガス低温成分の研究」五十嵐宣孝 (阪大) 第３章 : http://wwwxray.ess.sci.osaka-u.ac.jp/OskXrayTlabHP/home/papers/mthesis/FY2016/Igarashi_Mthesis_20170224.pdf
- 「X 線天文衛星 XMM-Newton を用いた大マゼラン雲内超新星残骸DEM L299の
  硬 X 線成分の起源の研究」
  埼玉大学杉本樹信氏第３章 : http://crystal2.heal.phy.saitama-u.ac.jp/OB/2013/sugimoto_shu-ron_final_2.pdf

データの検索

ESASky 2.2.1 (http://sky.esa.int)
- 絵を見ながらデータの検索が簡単にできる。
- 点源のチェック、アーカイブデータの有無、OBSIDなども取得できる。
- 視野、Chandra, 可視、電波の情報も一括で取得できる
ADS (http://adsabs.harvard.edu/abstract_service.html)
- 論文から天体を検索する王道。

ESASkyの使い方

天体を検索する。SS433などの天体名や、座標で検索する。
左上の Select Sky で、デフォルトで Optical で選択されているので、Soft X-ray で XMM-Newton EPIC color を選択する。
左上の渦巻きマーク (Image Observation) で、XMM-Newton を選択する。
視野内に含まれる観測情報、ObsIDなどが一覧で現れるので、所望のデータを選択する。

これで、XMM-Newton, Chandra, Suzaku, Integral などの基礎的な観測情報は取得できる。点源の情報なども、NEDやSimbadなどと連携していて、怪しげな点源はすぐに調べられる。

解析環境の構築

heasoft (https://heasarc.gsfc.nasa.gov/docs/software/lheasoft/)
- 概要
  - 標準的な解析ツール。これを入れずにX線解析はほぼ不可能。
  - ソースコードとバイナリの両方で配布している。32bit と 64 bit マシンのどちらでも動くバイナリを生成する仕様であるが、そのせいでコンパイルにコケることがある。
- linux へのインストール
  - ソースコードからコンパイルすれば、ハマることは少ないはず。ハマったら、make -n でコンパイルのコマンドを確認し、コンパイルできなかったライブラリを突き止めて、ライブラリの有無や変なオプションが指定されてないか等を確認する。xspec のローカルモデルを使う場合は、ソースコードからコンパイルしておく必要がある。
- mac へのインストール
  - mac へのインストールはハマる場合もそこそこある。OS毎に指示が違うので、マニュアルを丁寧に読もう。
  - https://heasarc.gsfc.nasa.gov/lheasoft/macosx.html
- anaconda python を使用時のコンパイルエラー
  - 手前味噌だが、これを参照 : https://qiita.com/yamadasuzaku/items/57829c607f6d5f15ce94
- perl のパスの設定
  - perl がデフォルトの #!/usr/bin/perl -w を使う設定であれば問題ないが、plenv を使っていると、emchain などの冒頭に、 #!/Users/syamada/.plenv/shims/perl -w を書かれていて、これだとなぜかシェルスクリプトから実行したときに、require: command not found のようなエラーがでて perl スクリプトが perl で実行されない状態に陥るので要注意。
sas (https://www.cosmos.esa.int/web/xmm-newton/sas-installation)
- 概要
  - XMM-Newton に特化した解析ツール。これがないとXMM-Newtonの解析はできない。
  - バイナリのみの配布。使ったことはないが仮想マシンでも配布している。
- Linux へのインストール
  - CentOS7 はOK、残念ながら debian には入らない。
- Mac へのインストール
  - El Captan では SIP (System Integrity Protection) を disable しなければならず、recovery mode で立ち上げて、csrutil disable をしないとダメ。
    - 詳細 : https://www.cosmos.esa.int/web/xmm-newton/sas-watchout-1600-sip
- 領域の指定
  - ds9 の領域ファイルを直接読み込んでくれるツールはない(もし違ったらご指摘ください)。例えば、evselect のオプションで、露わに書き下さないといけない。
    - evselect table='mos1_filt_time.fits' energycolumn='PI' withfilteredset=yes filteredset='mos1_filtered.fits' keepfilteroutput=yes filtertype='expression' expression='((X,Y) in CIRCLE('${cx}','${cy}','${rin}'))' withspectrumset=yes spectrumset='mos1_pi.fits' spectralbinsize=5 withspecranges=yes specchannelmin=0 specchannelmax=11999 > log/mos1_pi.log 2>&1
    - box だと ds9 と大きさが違う!?という話もあるようで、意図した領域が切り取られているかどうか確認しよう。
CCF (current calibration files)
- 較正ファイル一式。観測時期に対応した最新にキャリブレーションファイルを用いてデータをリプロセスして解析を行う。
- https://www.cosmos.esa.int/web/xmm-newton/current-calibration-files
- Valid CCF というのが、CCFの新しいファイルのみのサブセットで、600MB ある。通常はValid だけでよいと思われる。CCFは3GB。

データのダウンロード

例えば、obsid=0075140401 のデータを解析したい場合、下記の curl コマンドでデータをダウンロードする。

curl -o 0075140401.tar "http://nxsa.esac.esa.int/nxsa-sl/servlet/data-action-aio?obsno=0075140401"

上で、0075140401 の箇所を自分のOBSIDに変更すればよい。

ファイルの中には、

ODF - The ODF (raw) data files
PPS - The pipeline processed data products

が含まれており、PPSはパイプラインが解析してくれたプロダクトで、ざっとデータの概要を見たいときにはこの中のファイルを見ればよい。

普通は、生データ ODF をリプロセスして解析しなくてはならず、PPSは最初にちょっとだけ見る程度で、ほとんどの時間はODFと戯れることになる。

odf ディレクトリの中は、自動でファイルが解凍されてない場合があるので、tar で自力で中身を展開する。

tar -zxvf 0075140401.tar.gz 
tar -xvf 0881_0075140401.TAR

SASの初期設定

SASを使う上で、

https://www.cosmos.esa.int/web/xmm-newton/sas-thread-startup

は必ず読まないとダメ。SAS_DIR, SAS_PATHなどは自分がインストールしたSASの場所を示すように設定する。

例えば、heasoft と sas を /usr/local/astro 以下に入れた場合は、~/.bash_profile に下記の設定を書き込む。


export ASTRO=/usr/local/astro
export HEADAS=${ASTRO}/heasoft-6.22.1/x86_64-unknown-linux-gnu-libc2.17
alias heainit=". $HEADAS/headas-init.sh"
export PATH=${ASTRO}/ds9:$PATH
export PATH=${ASTRO}/xpa:$PATH
export PATH=${ASTRO}/wcstools-3.8.3/bin:$PATH

export SAS_DIR=${ASTRO}/sas/xmmsas_20160201_1833
export SAS_CCFPATH=${ASTRO}/sas/ccf
export SAS_ESAS_CALDB=${ASTRO}/sas/esas_caldb
export SAS_RAND_SEED=0
alias sasinit=" . ${SAS_DIR}/setsas.sh"

ややこしいが、CCF, ODF, CIF は重要な単語なので覚えておく。

Current Calibration Files (CCF) : キャルデータ
Observation Data Files (aka ODF) : 観測データ
Calibration Index File or CIF file : 使うべきキャルファイルを示したファイル

の３つは全部必要になる。

この時点で、sasinit を実行すると下記のようになる。

$ sasinit 

sasversion:- Executing (routine): sasversion  -w 1 -V 4
sasversion:- sasversion (sasversion-1.2)  [xmmsas_20160201_1833-15.0.0] started:  2018-07-18T11:28:08.000
sasversion:- XMM-Newton SAS release and build information:

SAS release: xmmsas_20160201_1833-15.0.0
Compiled on: Mon Feb  1 21:03:18 CET 2016
Compiled by: sasbuild@sasbld11
Platform   : RHEL6.6 64

SAS-related environment variables that are set:

SAS_DIR = /usr/local/astro/sas/xmmsas_20160201_1833
SAS_PATH = /usr/local/astro/sas/xmmsas_20160201_1833
SAS_CCFPATH = /usr/local/astro/sas/ccf

sasversion:- sasversion (sasversion-1.2)  [xmmsas_20160201_1833-15.0.0] ended:    2018-07-18T11:28:08.000

Do not forget to define SAS_CCFPATH, SAS_CCF and SAS_ODF

次に、SAS_CCF と SAS_ODF を定義する。ここでトリッキーな環境変数の付け替え操作が必要になるので、丁寧に説明する。

データをダウンロードした時点では、CIFファイル(キャルの必要情報が書かれたファイル)がない。これを最初に生成する。そのために、一旦、SAS_ODF がダウンロードしたデータのディレクトリを示すように設定する。

cd 0075140401/odf # odf ファイルのディレクトリに移動
export SAS_ODF=`pwd` # その場所にSAS_ODFを設定
cifbuild 
ls -thrl | tail -1 # 最新のファイルをチェック
> ccf.cif # ccf.cif ファイルが生成されたことを確認する。

次に、生成された ccf.cif を SAS_CCF が指し示すように設定し、odfingest を実行する。

export SAS_CCF=ccf.cif
odfingest 
ls -l | tail -1 # 最新のファイルをチェック
> 0881_0075140401_SCX00000SUM.SAS

SUM.SASファイルは、テキストファイルなので、中身を見て、コンテンツを確認できる。
これでSASの解析の準備が整った。

点源の解析

NASA/GSFCのABCガイドに実直にスクリプトを書くことにする。

https://heasarc.gsfc.nasa.gov/docs/xmm/abc/node8.html

この手順にそって、ほぼそのまま書きくだした格好にする。これに加えて、

https://heasarc.gsfc.nasa.gov/docs/xmm/abc/node15.html#Chap:epic-fit-xspec

に従い、スペクトルの groupting まで実行する。パイルアップの判定と、点源の判定は時間がかかるので割愛する。

パイルアップは、gradeの分岐比を使って、どの程度のパイルアップが起こっているかを判定するもので、あくまで判定ツールであって、補正するわけではないことに注意する。結構面倒くさいので、明るい天体でパイルアップを気にする人は解析が面倒になることを覚悟しよう (cf. [Suzaku pile up paper] (http://adsabs.harvard.edu/cgi-bin/bib_query?arXiv:1112.1844))

まずは、点源の中心座標を取得しよう。pps にイメージファイルがあるので、これを ds9 で見て、中心の detector 座標を取得する。

ds9 pps/P0075140401M2S002IMAGE_1000.FTZ # イメージファイルなので直接ds9で開ける

この例では、x,y = 28760, 35640 を中心座標をする。

次に、解析ディレクトルを作る。

mkdir -p ana/p1 # 解析ディレクトリの作成
cd ana/p1     
ln -s ../../0075140401/odf/*FIT . #イベントファイルのリンクを生成
ln -s ../../0075140401/odf/*SUM.SAS . # *SUM.SASファイルのリンクを生成
ln -s ../../0075140401/odf/ccf.cif .  # ccf.cif ファイルのリンクを生成

ここで、ABCガイドを実直に書き下した、

xmm_process_plotlc.sh 別リンクライトカーブをプロットするスクリプト
xmm_process_ql.sh 別リンクスペクトルまで作成する run スクリプト

を上記のディレクトリと同じ場所に置く。

スクリプトの中で、

export SAS_ODF=`/bin/ls -1 *SUM.SAS`
export SAS_CCF=ccf.cif

としているため、同じディレクトリに SUM.SAS, ccf.cif があることが前提となっている。

ここまで準備ができたら、あとは一気に解析を流すだけである。

スクリプトの引数は、cx cy rin rout の順番である。座標は physical で指定する。(ra, dec で指定したい人は、astropy などで変換しよう。)

中心座標(cx, cy) physical, (default: cx=25286.2, cy=18506.7)
点源を抽出する半径 = rin, physical (default=1200), 60 arcsec
bgdを抽出する半径 = rout, physical (default=3600), 180 arcsec


$ ./xmm_process_ql.sh 28760 35640 1200 3600 
cx       =  28760
cy       =  35640
rin      =  1200
rout     =  3600
[repro] process the MOS data
[repro] process the PN data
[everename] please check event file name
evem1 =  P0075140401M1S001MIEVLI0000.FIT
evem2 =  P0075140401M2S002MIEVLI0000.FIT
evepn =  P0075140401PNS006PIEVLI0000.FIT
[createDisplayImage] : mos1
[createDisplayImage] : mos2
[createDisplayImage] : pn
[stdfilt] : mos1
[stdfilt] : mos2
[stdfilt] : pn
[mklc] : mos1
[mklc] : mos2
[mklc] : pn
[mklc] : all
[mkgti] : create gtiset.fits
GTI is created using MOS1, if not, please modify here
[mkgti] : apply gtiset.fits into mos1 mos2 pn 
[srcDetect]
[srcDetect] : First, make the attitude file
[srcDetect] : create mos1-s
[srcDetect] : create mos1-h
[srcDetect] : create mos1-all
[srcDetect] : create mos1-all_bin_88
[srcDetect] : create mos1-all_bin_176
[srcDetect] : create mos2-s
[srcDetect] : create mos2-h
[srcDetect] : create mos2-all
[srcDetect] : create mos2-all_bin_88
[srcDetect] : create mos2-all_bin_176
[srcDetect] : create pn-s
[srcDetect] : create pn-h
[srcDetect] : create pn-all
[srcDetect] : create pn-all_bin_88
[srcDetect] : create pn-all_bin_176
[mkspec] : create src mos1
[mkspec] : create bgd mos1
[mkspec] : create mos1-all_bin_88 for src
[mkspec] : create mos1-all_bin_88 for bgd
[mkspec] : create src mos2
[mkspec] : create bgd mos2
[mkspec] : create mos2-all_bin_88 for src
[mkspec] : create mos2-all_bin_88 for bgd
[mkspec] : create src pn
[mkspec] : create bgd pn
[mkspec] : create pn-all_bin_88 for src
[mkspec] : create pn-all_bin_88 for bgd
[mkarea] : backscale for pi files
[mkrmfarf] : mos1
[mkrmfarf] : mos2
[mkrmfarf] : pn
[mkgrppha] : mos1
[mkgrppha] : mos2
[mkgrppha] : pn

とエラーがでずに流れればうまく走ったと思われる。

具体的な実行内容は、

./xmm_process_ql.sh


run()
{
setpass # or setpass $basedir $odffile
repro               # 6.1 Rerun the Pipeline
everename            # 6.1 (optional) re name event files
createDisplayImage   # 6.2 Create and Display an Image
stdfilt              # 6.3 Applying Standard Filters the Data
mklc                 # 6.4 Create and Display a Light Curve
mkgti                # 6.5 Applying Time Filters the Data
srcDetect            # 6.6 Source Detection with edetect_chain
mkspec               # 6.7 Extract the Source and Background Spectra
mkarea               # 6.10 Determine the Spectrum Extraction Areas
mkrmfarf             # 6.11 Create the Photon Redistribution Matrix (RMF) and Ancillary File (ARF)
mkgrppha              # 13. Fitting an EPIC Spectrum in XSPEC
}

のように、ABCガイドのセクション毎に関数化されて、これを順次実行する。

注意事項

点源の検出 (edetect_chain) は遅いからカットした。暗くてビンまとめが弱い場合は、10分くらいかかる。(どなたかコツご存知でしょうか？？あるいはそういうもんでしょうか。PSFで sliding してるからビンまとめが弱いと遅くなるんだと思ってるんですが。)
時間カットは、mos1 で rate < 6 を mos1 mos2 pn にも適用。

スペクトル解析

レスポンス、arf ファイルも生成されて、group min 20 でビンまとめされたスペクトルが生成されている。
プロットするには下記のように xspec で読み込む。

$ xspec

XSPEC12>data 1:1 mos1_pi_grp20.fits 
***Warning: Detected response matrix energy bin value = 0 (or neg).
     XSPEC will instead use small finite value (response file will not be altered).

1 spectrum  in use
 
Spectral Data File: mos1_pi_grp20.fits  Spectrum 1
Net count rate (cts/s) for Spectrum:1  5.978e-03 +/- 8.122e-04 (30.9 % total)
 Assigned to Data Group 1 and Plot Group 1
  Noticed Channels:  1-449
  Telescope: XMM Instrument: EMOS1  Channel Type: PI
  Exposure Time: 3.194e+04 sec
 Using fit statistic: chi
 Using test statistic: chi
 Using Background File                mos1_bgd_pi.fits
  Background Exposure Time: 3.193e+04 sec
 Using Response (RMF) File            mos1_rmf.fits for Source 1
 Using Auxiliary Response (ARF) File  mos1_arf.fits

XSPEC12>data 2:2 mos2_pi_grp20.fits 
***Warning: Detected response matrix energy bin value = 0 (or neg).
     XSPEC will instead use small finite value (response file will not be altered).

2 spectra  in use
 
Spectral Data File: mos2_pi_grp20.fits  Spectrum 2
Net count rate (cts/s) for Spectrum:2  4.733e-03 +/- 7.820e-04 (26.5 % total)
 Assigned to Data Group 2 and Plot Group 2
  Noticed Channels:  1-368
  Telescope: XMM Instrument: EMOS2  Channel Type: PI
  Exposure Time: 3.195e+04 sec
 Using fit statistic: chi
 Using test statistic: chi
 Using Background File                mos2_bgd_pi.fits
  Background Exposure Time: 3.195e+04 sec
 Using Response (RMF) File            mos2_rmf.fits for Source 1
 Using Auxiliary Response (ARF) File  mos2_arf.fits

XSPEC12>data 3:3 pn_pi_grp20.fits 

3 spectra  in use
 
Spectral Data File: pn_pi_grp20.fits  Spectrum 3
Net count rate (cts/s) for Spectrum:3  1.918e-02 +/- 1.909e-03 (21.3 % total)
 Assigned to Data Group 3 and Plot Group 3
  Noticed Channels:  1-379
  Telescope: XMM Instrument: EPN  Channel Type: PI
  Exposure Time: 2.712e+04 sec
 Using fit statistic: chi
 Using test statistic: chi
 Using Background File                pn_bgd_pi.fits
  Background Exposure Time: 2.712e+04 sec
 Using Response (RMF) File            pn_rmf.fits for Source 1
 Using Auxiliary Response (ARF) File  pn_arf.fits

XSPEC12>ignore 1-3:*-0.5 10.0-**

XSPEC12>setp e
XSPEC12>cpd /xs
XSPEC12>pl ld

これで次のように、MOS1(黒), MOS2(赤), PN(緑) の３つのスペクトルができることを確認しよう。

次に、

XSPEC12>setp com r x 0.1 10
1
XSPEC12>plot eff

として、有効面積をプロットする。

MOS1(黒)と MOS2(赤)はほぼ同じで、 PN(緑) はその２倍程度の有効面積になるはず。

xspec の使い方を知らない人は、下記をまずは一読しよう。

xspecの使い方
- 「すざく」first step マニュアルのスペクトル解析の章 : http://cosmic.riken.jp/suzaku/help/guide/fstep_web/node10.html
- 「すざく」衛星の日本語マニュアル。ds9, xspec, xronos の使い方は今でもほとんど変わっていない。
- ただし、「すざく」データの本格解析やリプロセスを行う場合は、必ず英語版の公式マニュアルを参照すること。2009年で日本語版のマニュアルの更新は終了し、英語版に一本化されたため。
  - http://cosmic.riken.jp/suzaku/help/guide/index_j.html
- Hongowiki
  - xspec12の使い方 : http://ytkyk.info/wiki/xspec/Xspec11ユーザのためのXspec12入門講座.html
  - 「すざく」のマニュアルは xspec11の頃なので、xspec12との違いはHongowikiで補完する
スペクトル解析の基礎事項
- 榎戸さんの資料
  - https://researchmap.jp/mu8wbt9s2-28502/?action=multidatabase_action_main_filedownload&download_flag=1&upload_id=5831&metadata_id=16414
  - 実験データの統計解析榎戸輝揚 (ver.2010-05-06@東大牧島研セミナー), フィットの基礎はこれで学ぶ。
- 星間吸収について
  - 海老沢先生の講義ノートの第５章ほか
    - http://www.isas.jaxa.jp/home/ebisawalab/ebisawa/TEACHING/2006UnivTokyo.html
    - http://www.isas.jaxa.jp/home/ebisawalab/ebisawa/TEACHING/20070122.pdf

イメージの確認

sas は ds9 形式の領域を受け付けない。また、長方形の場合は、ds9 とサイズがずれるという話もある。作成したら、必ず領域のカットがかかったイメージの ds9 で確認しよう。

領域フィルターのかかったソース領域は、
mos1_filtered_bin88.fits, mos2_filtered_bin88.fits ,pn_filtered_bin88.fits
を見て確認し、バックグランド領域は、
mos1_bgd_filtered_bin88.fits mos2_bgd_filtered_bin88.fits pn_bgd_filtered_bin88.fits
を見て確認しよう。bin88 は、デフォルトで適当に 88 程度のビンまとめがされるのためだけで、好みで変えてよい。

ライトカーブのチェック

xmm_process_ql.sh の中で、

mklc                 # 6.4 Create and Display a Light Curve

を走らせた時に、自動的に lcurve でライトカーブを生成して、mos1_filt_lc_mos2_filt_lc_pn_filt_lc.pdf として保存されている。

これで、所望の絵かどうか確認する。

論文用にエネルギーをカットしたライトーブや、時間ビンを変えたライトカーブを作成するしたい場合は、

https://www.cosmos.esa.int/web/xmm-newton/sas-thread-timing

に従って、ライトカーブを生成する。xmm_process_ql.sh の中で、

mklcps 500 2000 100 # 500eV-2000eV dt=100sec

という行がある。mklcps エネルギーの下限エネルギーの上限時間ビン、という使い方になる。

ライトカーブだけ作り直したい場合は、下記のように setpass だけ実行し、他はスキップして、mklcps だけ実行する。

run()
{
setpass # or setpass $basedir $odffile
#repro               # 6.1 Rerun the Pipeline
#everename            # 6.1 (optional) re name event files
#createDisplayImage   # 6.2 Create and Display an Image
#stdfilt              # 6.3 Applying Standard Filters the Data
#mklcps                 # 6.4 Create and Display a Light Curve
#mkgti                # 6.5 Applying Time Filters the Data
#srcDetect            # 6.6 Source Detection with edetect_chain
#mkspec               # 6.7 Extract the Source and Background Spectra
#mkarea               # 6.10 Determine the Spectrum Extraction Areas
#mkrmfarf             # 6.11 Create the Photon Redistribution Matrix (RMF) and Ancillary File (ARF)
#mkgrppha              # 13. Fitting an EPIC Spectrum in XSPEC

# extra process 
mklcps 500 2000 100 # 500eV-2000eV dt=100sec, sas-thread-timing : EXTRACTION AND CORRECTION OF AN X-RAY LIGHT CURVE FOR A POINT-LIKE SOURCE

}

デフォルトで 100秒の時間ビンで、0.5keV - 2 keV のライトカーブを生成する。
作成したライトカーブは、mos1_net_filtered_lc_500_2000_100_mos2_net_filtered_lc_500_2000_100_pn_net_filtered_lc_500_2000_100.pdf に保存されるので、結果を確認してほしい。

ソースとバックグランドの領域を微妙に変えたい場合

下記、４通りの初期値の与え方がデフォルトで用意してあるので、適宜修正してお使いください。

./xmm_process_ql.sh cx cy rin rout
./xmm_process_ql.sh cx cy rin rout rin2 (when rin2 is for background)
./xmm_process_ql.sh cx cy rin rout rin2 (when rin2, cx2, cy2 is for background)
./xmm_process_ql.sh cx cy rin rout rin2 binsize (when rin2, cx2, cy2 is for background, binsize is group min value)

その他の方法

xmmextractor
- https://www.cosmos.esa.int/web/xmm-newton/sas-thread-xmmextractor
- ra, dec を XML で指定すれば、点源の解析を一挙にやってくれるツール。
- 領域の半径、バックグランド領域の指定方法が不明なのが論文化の時に問題になりそう。(わかる方いたら教えてください。)

真面目に解析する前に知っておくこと。

バックグラウンド源

スペクトル解析でもイメージ解析をする上でも、バックグラウンドにどんな種類があるのかは把握しておこう。

バックグラウンド源
- ~ 数100keV程度の陽子
- 宇宙線陽子
- 電気的ノイズ
- 宇宙背景放射
- 軟X線放射
  - Local bubble
  - Galactic disk
  - Galactic Halo
  - Solar Wind Charge Exchange (SWCX)
  - 視野外からの一回反射成分
  - CCDの時間外成分

[参考ページ] : https://www.cosmos.esa.int/web/xmm-newton/epic-background-components

広がった天体の解析

イメージ解析

広がった天体の解析では、エネルギーバンドごとに分けたイメージを作成し、データをよく見ることから解析が始まる。ここでは、SASのスレッド

「CREATION OF EPIC BACKGROUND SUBTRACTED, EXPOSURE CORRECTED IMAGES」
- https://www.cosmos.esa.int/web/xmm-newton/sas-thread-esasimage

に従ってイメージを生成する。このスレッドの通りに実行するスクリプトが

xmm_process_plotlc.sh 別リンク: バックグランドを引いて、Exposure補正したイメージを生成するスクリプト (一戸悠人君によって作れらました。)

である。ただし、cosmic ray flux はフィットが必要なので、自動化には含まれてにない。気になるレベルの解析をする場合は、10. Fit the spectral data to determine the soft proton contamination parameters. で soft proton をちゃっと決める部分もやろう。

実行するには、まず作業ディレクトリを作成し、

mkdir -p image
cd image

イベントファイル、SUM.SASファイル、ccf.cifファイルのリンクを取得し、

ln -s ../../0075140401/odf/*FIT .
ln -s ../../0075140401/odf/*SUM.SAS .
ln -s ../../0075140401/odf/ccf.cif .

xmm_process_image_each.sh を実行する。

./xmm_process_image_each.sh

エラーがでずに終了していたら、バックグランドを引かれて、exposure を補正された３つのイメージが出力されている。

0.4-0.75 keV : adapt-400-750-all.fits
0.75-1.3 keV : adapt-750-1300-all.fits
2.0-7.2 keV : adapt-2000-7200-all.fits

ds9 で生成されたイメージを確認する。

ds9 -rgb -red adapt-400-750-all.fits -green adapt-750-1300-all.fits -blue adapt-2000-7200-all.fits

そうすると、このような３色のイメージ生成される。

WCS空間での画像の足し方

WCS空間で画像を足したり、exposure 補正する方法は下記を参照。
https://qiita.com/yamadasuzaku/items/62f55503c29d6a6cc1aa

OBSIDのみで自動イメージ解析するスクリプト

下記にOBSIDだけ入れると、上記のイメージ解析を一挙に行うスクリプトを置いた。
http://www-x.phys.se.tmu.ac.jp/%7Esyamada/syamadatools/xmm-newton/xmm_process_autorun.sh

./xmm_process_autorun.sh $obsid

で実行すると、データのダウンロード、展開、イメージ解析まで一挙にやってくれる。

radial profile の作り方

点源がどうか、広がっているかどうかを radial profile　から確認する方法を示す。

xmm_process_radialprofile.py : SKY座標を指定して、その座標点まわりの radial profile を作成する。
anaconda2系でastropyなどがインストールされていることが前提。

# ファイルのリストを生成
ls adapt-750-1300-all.fits > xmm_adapt.list 
# SKYX=218, SKYY=152 半径=10 で実行する。詳細は -h 
./xmm_process_radialprofile.py xmm_adapt.list -x 218 -y 152 -s 10

成功すると下記のような絵が、fig_image_xy_adapt-750-1300-all というディレクトリの下に生成されている。

fitsのイメージをROOTのTH2Dに変換

ROOTユーザーの場合、fits ファイルをROOTのTH2Dに変換して解析したいこともある。
そのためのスクリプトも用意しておいた。

fitsimage2root.py : FITSのイメージをROOTのTH2Dに変換する。
- pyROOTがインストールされていることが前提。
- オプションで、指定したSKY領域のイメージのみをダンプすることも可能。イメージが大きすぎる場合にカットできて便利。詳細は -h オプションを参考。

# ファイルのリストを生成
ls adapt-750-1300-all.fits > xmm_adapt.list 
# オプション無しで実行すると全イメージを変換する。詳細は -h 
./fitsimage2root.py xmm_adapt.list

生成されたROOTファイルはTBrowser で確認しよう。使い方は、

%root -l adapt-750-1300-all.root                   
root [0] TBrowser b

で TBrowser を立ち上げて、オプションを"COLZ"などを指定して描画する。

選択した領域のヒストグラムなどROOTの機能を使えば簡単にできるので、ROOTユーザーには便利かもしれない。