RHESSI データの解析方法 II

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

(0) 観測データの取得

ディレクトリー構造は YYYY -> MM -> DD となっており、 欲しい日付までたどり着いたら、hsi_yyyymmdd_hhmmss_vvv.fits というfileがあるので、欲しい時間帯のfileをダウンロードする。
＊観測fileの形式はFITSです。
＊将来的にはフレアリストに自動でアクセスして、データ取得の手間は なくなる予定です。

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

IDL> obj_ltc=hsi_lightcurve()
obj_ltcというobjectを定義。

IDL> obj_ltc->plot, obs_time_interval=['2002/04/04 15:00:00','2002/04/04 16:00:00']
ライトカーブを描かせたい時間帯を指定。更に、エネルギー帯を指定したりするには、
IDL> obj_ltc->plotman, obs_time_interval=['2002/04/04 15:00:00','2002/04/04 16:00:00'], $
ltc_energy_band=[25,50,100,150],ltc_time_resolution=1.0
ここでplotではなく、plotmanにしたのはplot managerを立ち上げるため (色別で描いてくれたりする)。
ltc_time_resolutionはplotする時間を1秒間隔にする。

(2) 像合成

(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 を参照してください。