Xspecで自分のモデルをテーブルとして入れる

このエントリーをはてなブックマークに追加

XspecでX線データの解析をしていると、自分で作成したモデルで天体スペクトルをフィットしてみたくなることがあります。 とてもありがたいことに(開発メンバーに感謝しましょう) Xspec version 12からは、かなり簡単に、自分のモデルを関数として 入力するための枠組みが整っています。このページでは、実際にそれを利用してモデルを構築する方法を紹介します。

1. 関連するページ

独自モデルに関連するページは以下のようなものがあります。

  1. Xspecのマニュアルの解説ページが大変参考になります。
  2. 関数モデルではなく、スペクトル情報をFITSテーブルとして与える方法については、 Xspecで自分のモデルをテーブルとして入れるを参照してください。
  3. 山田くんによる、相対論を考慮したKYモデルの使い方のページでは、 他の研究者が作成した独自モデルをダウンロードしてきて組み込む方法を 解説してくれています。

2. 利用するコマンドと作業の流れ

関数モデルを入れる際には、initpackageとlmodというコマンドを利用します。 initpackageコマンドは、コマンドラインからも利用できますが、xspec内から 呼ぶこともできます。

作業の流れは以下の通りです。 0. '''ソースコード版のHeasoftをインストール(バイナリ版ではなくて)''' 自分の関数をXspecが読み込めるダイナミックリンクライブラリとしてコンパイルするときに、 Xspecのパッケージに含まれるヘッダファイルが必要なので。 1. '''新規モデルの説明ファイルを作成''' 追加しようとしているモデルの名前やパラメータ数、パラメータの範囲などを指定した テキストファイルを、指定のフォーマットで作成します。ファイル名は何でもよいのですが、 よくlmodel.datというのが用いられるようです(eqpairの場合など)。 2. '''関数モデルの準備''' 追加する関数モデルのコードを記述します。 計算すべきエネルギー範囲や、各エネルギービンのステップ、モデルパラメータなどが 引数として与えられるので、それに対応したモデルフラックスを計算して、 配列(もしくはC++のvalarray)として返すような関数をつくることになります。 Fortran、C/C++で記述することができます。おすすめはC++。 3. '''共有ライブラリとしてビルド(initpackage)''' 作成した関数を、Xspecがプログラム内から呼び出せるように、共有ライブラリとしてコンパイルします。 Xspecではこうして作成される共有ライブラリのことを「パッケージ」と呼ぶようです。 上で作成したlmodel.datと関数がただしければ、initpackageコマンドが、Makefileの作成、ビルドまで自動的に行ってくれます。たとえば、今回作成したモデル名が「mymodel」で、パッケージ名 として「mypackage」を指定した場合は、libmypackage.soという共有ライブラリが作成されます。 4. '''Xspecでの読み込み''' 作成したモデルをlmodコマンドで読み込みます。上の例で、libmypackage.soという パッケージ(共有ライブラリ)が作成されている場合は、lmod packageとすればOKです。 パッケージをサーチするフォルダパスは、~/.xspec/Xspec.initの中で指定しておく 必要があります。パッケージの読み込みが完了すると、modelコマンドでモデルをしていするときに、 追加したモデルが利用可能になります。上の例では、

    model mymodel+pow

みたいな記述が可能になるわけです。

3. 具体的な説明

以下、まず作業上の注意点を述べてから、具体的な手順を解説します。

3.1 注意すべきところ

Unable to find a resource file appropriate for HEADAS component Xspec

なぜか、initpackage(の内部で利用されているhmake)は、使用しているファイルのパスに シンボリックリンクがあるとうまく動作しないことがあります。とくに、HEADAS環境変数 (シェルで、echo $HEADASとしたときに表示されるパス)にシンボリックリンクがあると、 initpackageでのビルドの際に、

hmake: Unable to find a resource file appropriate for HEADAS component Xspec
hmake: in the following path: /users/yuasa/tmp/ippsr/sources:/users/yuasa/tmp/ippsr/sources/BUILD_DIR:
/usr/local/headas/6.8/headas/linux/bin:/usr/local/headas/6.8/headas/Xspec/linux/BUILD_DIR

みたいなエラーがでます。上の例では、エラーの原因はHEADAS環境変数を指定するときに

export HEADAS="/usr/local/headas/6.8/headas/x86_64-unknown-linux-gnu-libc2.3.4"

と指定すべきところを、x86_64-unknown-linux-gnu-libc2.3.4が長いので、linuxというシンボリックリンクを 作成して、

export HEADAS="/usr/local/headas/6.8/headas/linux"

としていたことが原因でした。/usr/local/headas/6.8/headas/x86_64-unknown-linux-gnu-libc2.3.4を 指定するようにすると、上記のエラーはでなくなりました。

3.2 新規モデルの説明ファイルを作成

以下のようなフォーマットで、説明ファイルを作成します。 ファイル名は何でもよいのですが、lmodel.dat(local modelの意)がよく使われるようです。

1行目でモデル全体の概要を記述し、2行目以降は、各パラメータについて指定を行います。

(model name) (number of model parameters) (model lower limit E) (model upper limit E) (model type) (C/C++ function name) (flag for model variance) (flag for recalculation)
(parameter1 name) (unit) (default value) (hard min) (soft min) (soft max) (hard max) (fit delta)
(parameter2 name) (unit) (default value) (hard min) (soft min) (soft max) (hard max) (fit delta)
...

モデルタイプは、add, mul, mix, or con, acnのどれかになります。 1行目の最後のふたつのフラグの説明は、オフィシャルの説明を参照してください。 通常は「0 0」と記述すれば大丈夫でしょう。 パラメータの単位のところは、ユーザの便宜のためのに表示されるものなので、 わかりやすいものをいれればOKで、無次元の場合は、" "を記述します。 各リミット(hard minなど)やfitting deltaについては、newparコマンドのヘルプをみてください。

具体的には、以下のようになります。これ以降、モデル名は、ippsrとします(Intermediate Polarの Post Shock Regionからの放射モデルという意味)。この例では、ふたつのパラメータを定義しています。 関数名はC_ippsrとしました(冒頭のC_は、XSpec的にいうと、C++で記述している、という意味ですが、 あってもなくても大丈夫です)。 モデルタイプは、additiveモデル(power lawやgaussianなどと同じ、そのコンポーネント自信がフラックスを 表現するモデル)になります。Massは白色矮星の質量を太陽質量単位で、ZFeは降着ガスのアバンダンスを 意味します。

ippsr 2 0 1e5 C_ippsr add 0 1
Mass Msun 0.6 0.1 0.1 1.39 1.39 0.01
ZFe " " 0.5 0.1 0.1 1.5 1.5 0.01

3.3 関数モデルの準備

lmodel.datで記述した関数を実装します。 ファイル名は、C_(モデル名).ccとして、保存場所はlmodel.datと同じにします。 prefixとして付与するC_とかc_が、引数の型を決めます(大文字のCだと、C++的な引数になり、小文字のcだとC言語的な引数になります)。

引数の型などについては、オフィシャルの説明を参照してください。

ここでは、C++のコードを例に説明します。

#include <xsTypes.h>
#include <iostream>
using namespace std;

extern "C" {
void ippsr(const RealArray& energy, const RealArray& parameter, int spectrum, RealArray& flux, RealArray& fluxError, const string& init){
    flux.resize(energy.size());
    //実際の計算コード

}
}

まず、xsType.hのincludeが必要です。この中で、Real型とRealArray型が定義されています。 Realはdoubleのようなもの、RealArrayはSTLのvalarrayのようなものと思ってください。 計算すべきenergyのビン情報が、energyに入っていますので、それに対応するフラックスを fluxに入れて返します。energyの単位はkeVです。 モデルのバリアンスを入れたい場合はfluxErrorにそれらを入れておきます。 initは、ユーザから実行時にパラメータを渡す機能(たとえば計算に使用するデータファイルを xspec内から指定する)などに利用できますが、必要ない場合は、宣言だけしておいて使用しなくてOKです。 実際の計算コードのところで、何らかの計算をして、fluxにフォトンフラックスを つめていってください。

重要なポイント1 : fluxアレイの単位

fluxアレイに入れるべきものは、オフィシャルの説明に あるように、

The output flux array for an additive model should be in terms of photons 
cm2 s-1 (not photons cm2 s-1 keV-1) i.e. it is the model spectrum integrated over the energy bin.

です。日本語で繰り返すと、fluxアレイに入る数値は、そのビン幅におけるフォトンフラックスをビン幅について積分したもの(photon cm^-2 s^-1)であり、ビン幅で規格化していないもの(つまり、ビン幅 keVで割らない)です。誤ってビン幅で割ってしまうと、ぜんぜん違った結果になってしまうので注意。自分で入れたモデルが正しいかどうか、xspec内で、plot modelして確認しましょう。

重要なポイント2 : fluxアレイのサイズ

additiveモデル(現在説明しているモデルタイプ)では、fluxとfluxErrorのサイズは、 関数呼び出し時は0になっている事に注意してください。実際の計算を行う前に、 上記の例のように、resizeメソッドを使って、energyアレイと同じサイズに設定する 必要があります(ちょっと面倒で、忘れると実行時にsegmentation faultになる)。 multiplicativeモデルでは、計算にしようするスペクトルデータが入っているので、 resizeの必要はありません(元からenergyアレイと同じサイズなので、resizeしても 何もおきません)。

4. 具体的なサンプル

powerlawと同じ関数系を入れる場合のサンプル(lmodel.datとc_testmodel.ccファイル)を含んだzipパッケージを載せておきます。

Download: Media:testmodel.zip

アーカイブを展開して、

cd testmodel
initpackage testmodelpackage lmodel.dat .
hmake

とすると、libtestmodelpackage.soが生成されます。

エディタでXspecの初期設定ファイルを開き、 nano \$HOME/.xspec/Xspec.init に LOCAL_MODEL_DIRECTORY: /home/yuasa/tmp/testmodel 等、libtestmodelpackage.soが入っているフォルダを指定します。

以下のように、Xspecを起動してlmod testmodelpackageコマンドでロードすると、model testmodelという形で自分で作ったモデルが使用可能になります。

        XSPEC version: 12.8.1
    Build Date/Time: Fri Aug 30 18:50:01 2013

 Cannot find ATOMDB_VERSION in Xspec.init,  using default : (2.0.2)
XSPEC12>lmod testmodelpackage
Model package testmodelpackage successfully loaded.
XSPEC12>model testmodel

Input parameter value, delta, min, bot, top, and max values for ...
            2.1       0.01(     0.021)        0.1        0.1          5          5
1:testmodel:index>
              1       0.01(      0.01)          0          0      1e+20      1e+24
2:testmodel:norm>

========================================================================
Model testmodel<1> Source No.: 1   Active/Off
Model Model Component  Parameter  Unit     Value
 par  comp
   1    1   testmodel  index               2.10000      +/-  0.0          
   2    1   testmodel  norm                1.00000      +/-  0.0          
________________________________________________________________________

testmodel.zip内のそれぞれのファイルの中身は次の通りです。

lmodel.dat

testmodel 1 0 1e5 c_testmodel add 0 1
index " " 2.1 0.1 0.1 5 5 0.01

c_testmodel.cc

#include <vector>
#include <xsTypes.h>
#include <iostream>
#include <cmath>
using namespace std;

extern "C" {
    void testmodel(
        const Real* energy, //0--Nflux+1, in keV
        int Nflux, 
        const Real* parameter, 
        int spectrum, 
        Real* flux, 
        Real* fluxError, 
        const char* init
    ){
        double Gamma=parameter[0];
        double B=1/(-Gamma+1);

        using namespace std;
        /* cout << "Energy[0]= " << energy[0] << endl;
         *  cout << "Gamma    = " << Gamma << endl;
         */

        for(size_t i=0;i<Nflux;i++){
            double eLow=energy[i];
            double eUp=energy[i+1];
            flux[i]=B*( pow((eUp),-Gamma+1) - pow((eLow),-Gamma+1) );
        }

    }

}

4.1 共有ライブラリとしてビルド(initpackage)

作成した関数モデルをビルドします。 lmodel.datがあるフォルダ内でXspecを起動して、以下のコマンドを入力します。

XSPEC12>initpackage mypackage lmodel.dat .

第一引数は、パッケージ名で、共有ライブラリはlib(パッケージ名).so/dyldとなります。 第二引数に、lmodel.datを指定します。 第三引数に、作業フォルダ(上の例では、lmodl.datと同じフォルダ)を指定します。 自動的にMakefileが作成され、hmakeを用いたビルドが実行されます。 作業が完了すると、libmypackage.soが作成されているはずです。

4.2 Xspecでの読み込み

まず、~/.xspec/Xspec.initをテキストエディタで開いて、 LOCAL_MODEL_DIRECTORY:という変数に、上の作業で 共有ライブラリを作成したフォルダパスを指定しておきます。

#エディタで開く
nano ~/.xspec/Xspec.init

LOCAL_MODEL_DIRECTORY: (共有ライブラリを含むフォルダへのパス)

例:
LOCAL_MODEL_DIRECTORY: /users/yuasa/ip/models/ippsr_2010022

この設定後、作成した共有ライブラリを以下のようにXspecで読み込み、 以下のコマンドを実行して、

XSPEC12>lmod mypackage
Model package mypackage successfully loaded.

となれば、そのパッケージに含まれる関数モデルが利用できるようになります。 この場合は、ippsrというモデル名だったので、以下のようにすると利用できるわけです。

XSPEC12>model ippsr

Input parameter value, delta, min, bot, top, and max values for ...
            0.6       0.01        0.1        0.1       1.39       1.39
1:ippsr:Mass>0.8
            0.5       0.01        0.1        0.1        1.5        1.5
2:ippsr:ZFe>0.6
              1       0.01          0          0      1e+24      1e+24
3:ippsr:norm>1
energy bin size : 51

========================================================================
Model ippsr<1> Source No.: 1   Active/Off
Model Model Component  Parameter  Unit     Value
 par  comp
   1    1   ippsr      Mass       Msun     0.800000     +/-  0.0          
   2    1   ippsr      ZFe                 0.600000     +/-  0.0          
   3    1   ippsr      norm                1.00000      +/-  0.0          
________________________________________________________________________

この状態で、

XSPEC12>cpd /xs
XSPEC12>plot model

とすると、モデルスペクトルが(検出器レスポンスがかかる前の形で)表示されるので、 自分が意図したものになっているか(スペクトル形とnormalization)を確認してください。

パラメータを入力した直後にSegmentation Faultが発生する場合は、 関数モデル内部で不正なメモリアクセスを行っている事が問題なので ソースコードを再検討してください。

5. HASOFT 6.14でのsymbol not found errorへの対処

2013-10-22時点での記述

Mac上のHeasoft 6.14の環境でinitpackageしローカルモデルをビルドした場合、lpack_(module名).cxxにUnload()とSafeUnload() という、Tclで必要な関数が見つからないといってエラーになりました。

XSPEC12>load /Users/yuasa/work/ip/model_calculation/model/
model_20100726/construct_ippsr2/libippsrpackage.dylib

dlsym(0x7decdf00, Ippsrpackage_Unload): 
symbol not founddlsym(0x7decdf00, Ippsrpackage_SafeUnload): symbol not found

とりあえず、lpack_ippsrpackage.cxxの冒頭に、

extern "C" int Ippsrpackage_SafeUnload(Tcl_Interp* tclInterp){}
extern "C" int Ippsrpackage_Unload(Tcl_Interp* tclInterp){}

という空の関数のエントリを作成してから、再度

hmake

すると、Xspec内でloadできるようになるので、当面これでしのいでいます。