RHESSI データの解析方法

このマニュアルは、 RHESSI解析マニュアルに準じて書きました。
Object Syntaxで作られているので、そのように書くようにしました。

(0) 観測データの取得

http://hesperia.gsfc.nasa.gov/hessidata/
ftp://hercules.ethz.ch/pub/hessi/data/
でも取得できる。

ディレクトリー構造は YYYY -> MM -> DD となっており、 欲しい日付までたどり着いたら、hsi_yyyymmdd_hhmmss_vvv.fits というfileがあるので、欲しい時間帯のfileをコピー(ダウンロード)する。
＊観測fileの形式はFITS。

(1) ライトカーブを描く

パラメータを指定
IDL> obj_ltc -> set,obs_time_interval=['2002/07/23 00:00:00','2002/07/23 01:10:00']
*この時間帯のデータを読み出す(少し長めに)
IDL> obj_ltc -> set,time_range=['2002/07/23 00:00:00','2002/07/23 01:10:00']
*ライトカーブを描かせたい時間帯
IDL> obj_ltc -> set,ltc_energy_band=[25,50,100]
*エネルギー帯を指定(この場合、25〜50、50〜100keVの2バンドを指定) IDL> obj_ltc -> set,ltc_time_resolution=1.0
*時間分解能(Bin)を指定(この場合1秒)

いよいよプロット!
IDL> obj_ltc -> plot
*ひたすら待つ!
また、ここで代わりに
IDL> obj_ltc -> plotman
とすると、plot managerが立ち上がる(色別で描いてくれたりする)。

これらの命令を、
IDL> obj_ltc -> plotman, obs_time_interval=['2002/07/23 00:00:00','2002/07/23 01:10:00'], $
IDL> ltc_energy_band=[25,50,100],ltc_time_resolution=1.0
と一気に書いてもOK

(2) 像合成

パラメータを指定
IDL> obj_im -> set,obs_time_interval=['2002/07/23 00:00:00','2002/07/23 01:10:00']
*この時間帯のデータを読み出す(少し長めに)
IDL> obj_im -> set,as_roll_solution='PMT',ras_time_extension=[-500,500]
IDL> obj_im -> set,energy_band=[25.,40.]
*エネルギー帯を指定(この場合、25〜40keV)
IDL> obj_im -> set,image_dim=[64,64]
*出来上がりの画像のサイズ(ピクセル数)を指定
IDL> obj_im -> set,pixel_size=[2,2]
*出来上がりの画像のピクセルサイズ(arcsec)を指定
IDL> obj_im -> set,xyoffset=[-900,-235]
*画像の中心座標、太陽中心を(0,0)として、arcsec単位で
IDL> obj_im -> set,time_range=['2002/07/23 00:20:00.000','2002/07/23 00:25:00.000']
*像合成する時間域を指定
IDL> obj_im -> set,det_index_mask=[0,0,0,1,1,1,1,1,0]
*使用するグリッドの番号。順に、1番、2番、、、9番。また、0はそのグリッドを使用しない、という意味
IDL> obj_im -> set,as_spin_period=4.0
IDL> obj_im -> set,image_algorithm='clean',clean_niter=100
*像合成プログラムを指定

これらの命令を、
IDL> obj_im -> set,energy_band=[25.,40.],image_dim=[64,64], $
IDL> pixel_size=[2,2],xyoffset=[-900,-235], $
IDL> time_range=['2002/07/23 00:21:00.000','2002/07/23 00:23:00.000'], $
IDL> det_index_mask=[0,0,0,1,1,1,1,1,0],as_spin_period=4.0
IDL> obj_im -> set,image_algorithm='clean',clean_niter=100
と一気に書いてもOK

いよいよ像合成!
IDL> obj_im -> plot
IDL> obj_im -> plotman
*合成した2次元画像を表示
IDL> im=obj_im -> getdata()
*合成した2次元画像を変数 im に入れる
IDL> map_hsi = obj_im -> make_map()
*map に変換

(3) 像合成のプログラム

1. back projection ; まずはざっと放射源をチェックするとき
   (imageのpeakの場所をpoint sourceとして全photon数をばらまく)
   IDL> hessi_image,pars,im,/plot
   するとこんな感じ

2. clean (back projection のcleaning)
   IDL> hessi_image,im,pars,/plot,image_algorithm='clean'
   するとこんな感じ
   

   使用するコリメータを変え、細かい構造まで見てみるには、(時間かかる)
   IDL> hessi_image,im,pars,/plot,image_algorithm='clean', det_index_mask=[0,0,1,1,1,1,1,1,1]
   するとこんな感じ

3. MEM SATO (佐藤さん作成のMEM) ; sharp だけど時間かかります
   なんだかまだうまくいっていないみたい
   (Sato et al. 1999, PASJ 参)
   MEMって？ … Maximum Entropy Method
   IDL> hessi_image,im,pars,/plot,image_algorithm='mem sato'
   するとこんな感じ
   

   (ここで補足)
   小さいピッチのコリメータを使った方が小さい構造まで追えます。 しかし、MEM SATOを計算するにはとっても時間がかかるので、 デフォルトでは番号にして4から8のコリメータ (0〜8の、全部で9つのコリメータがあり、番号が小さい方がピッチが小さい) を使うようにしてあります。
   使用するコリメータを指定するときは、
   IDL> hessi_image,im,pars,/plot,image_algorithm='mem sato', det_index_mask=[0,0,0,0,1,1,1,1,1]
   (この例では4,5,6,7,8のコリメータを使用する)
   

4. MEM VIS (MEM SATOより時間かかる)
   まだうまくいっていない
   (light curveではなく、visibility dataでMEMを使用している)
   IDL> hessi_image,im,pars,/plot,image_algorithm='mem vis'
   するとこんな感じ
   MEM SATOのようにdet_index_maskも使えます。
   

5. PIXON (更に時間かかる ; マニュアルによると750MHZ Pentium IIIのPCで 3時間半かかったらしい、私はもっとかかったみたい)
   まだうまくいっていない
   (Puetter and Pina, Experimental Astron., 3, 293; Pina and Puetter, 1993, PASP, Metcalf et al., 1996, ApJ 参)
   IDL> hessi_image,im,pars,/plot,image_algorithm='pixon'
   するとこんな感じ

6. Forward Fitting (結構早くてロバスト)
   たくさんの2D Gaussianの重ね合わせ
   IDL> hessi_image,im,pars,/plot,image_algorithm='forwardfit'
   これはGaussian 1つで計算。
   IDL> hessi_image,im,pars,/plot,image_algorithm='forwardfit', ff_n_gaussians=2
   でGaussianの数を指定する(この場合は2つ)。
   するとこんな感じ
   

ここには一例を挙げました。 詳しくはhttp://hesperia.gsfc.nasa.gov/ssw/hessi/doc/hsi_params_image_standard.htm を。

• image_algorithm ; 像合成時のプログラムを指定。(clean、MEM **等)

• det_index_mask ; 9x3の行列。1〜9はコリメータの番号、1〜3はハーモニクスの 番号(0;fundamental,..2;the third harmonics)。

• energy_band ; エネルギーバンド。デフォルトは[12.,25.]。
  例)
  IDL> hessi_image,im,pars,/plot,image_algorithm='clean',energy_band=[30.,80.]
  この例では[30 keV〜80 keV]のエネルギーバンド。 するとこんな感じ

• image_dim ; 像合成時のピクセル数、デフォルトは[64x64]。
• pixel_size ; 像合成時のピクセルの大きさ、デフォルトは[1.,1.](arcsec)。
  例)
  IDL> hessi_image,im,pars,/plot,image_algorithm='clean',image_dim=[32,64]
  するとこんな感じ
  IDL> hessi_image,im,pars,/plot,image_algorithm='clean',pixes_size=[2.,1.]
  するとこんな感じ

• xyoffset ; 合成する像の中心の座標を指定。
  太陽中心から、[600.,200.]のようにarcsecで入れる。

• obs_time_interval ; 観測があるかないかに関わらず、 とにかく興味ある時間帯を指定。
• time_range ; 像合成する時間帯。
  どちらも['yyyy/mm/dd hh:mm:ss','yyyy/mm/dd hh:mm:ss']で指定。

• image control parametersを返したいときは、
  IDL> hessi_image,im,pars,/plot
  で、parsにパラメタの情報が入り、それを読むには、
  IDL> help,pars,/st
  とすれば良い。
  

• RHESSIデータの解析に便利なコマンド集を(別のwindowで)呼び出す。
  IDL> hessi_list
  

• RHESSI解析用softwareは、release areaとdevelopment areaに分かれています。
  release area : test済みのsoftware
  development area : 現在もtest中のsoftware
  RHESSIデータ解析用softwareは現在も発展中であり、今はdevelopment areaに PATHが通っています(将来的にはすべてrelease areaにsoftwareが置かれ、 PATHもそちらに通るようになります)。これらのareaの切替えは、
  IDL> hessi_version,/dev
  IDL> hessi_version,/release
  で行ないます。
  

(4) スペクトル解析

(指定する必要のあるキーワード達)

• binning
  IDL> hsi_spectroscopy_list
  スペクトル解析の方法をリストアップしてくれる。
  例えば、6番のスキームを使うなら、
  IDL> sp -> set,sp_energy_binning=6
  

• 時間間隔(積分時間)
  IDL> sp -> set, sp_time_interval=10
  とすれば、10秒間隔でたくさんのprofileを取得出来る。
  IDL> sp -> set, sp_time_interval=['2002/04/04 15:36:10','2002/04/04 15:36:14']
  とすれは、積分する時間を指定できる。 例えば、ピークの時間だけスペクトルを取りたい時などに使う。

• detector
  現在は各コリメータのfundamental modulationに対してしか計算できないので、 そのように指定します。
  IDL> sp -> set, det_index_mask=bytarr(9)+1b
  

(さらに補助的だけど大事なキーワード)

• IDL> spec = sp -> getdata(sum_flag=0)
  コリメータ毎のスペクトルを得る (これを指定しないと、全コリメータの値の総和になる)。

• IDL> spec = sp -> getdata(sp_data_unit='Flux')
  フラックスにする。 [phorons /s /cm^2 /keV]
  
• IDL> spec = sp -> getdata(sp_data_unit='Counts')
  カウントにする(デフォルトに戻す)。
  
• IDL> spec = sp -> getdata(sp_data_unit='Rate')
  [phorons /s /keV]

• IDL> spec = sp -> getdata(/sp_data_structure)
  スペクトルに関するデータをstructureで返す (プロットするときに便利なのでやっておいた方が良い)。
  
• IDL> spec = sp -> getdata(/sp_semi_calibrated)
  これは、grid transmission、attenuator state、detector response等々 といったものに対する補正のためなので、 必ず行なってください！

＊GUIやSPEXでスペクトル解析をするなら http://hesperia.gsfc.nasa.gov/~dennis/spectroscopy/first_steps.htm を参照してください。