---
lang:ja-jp
###### tags: `open`
---

# NuSTARの解析（藤原メモver）
++[X線スペクトル解析](/wOsMOU-3QCC6rkHt3kdfuA)++
++[衛星ごとの処理(まとめ)](https://hackmd.io/@ystl/rkOK0QvPS)++

++[Xspec($\in$ heasoft)の導入](/5sw6ekv1RwSyckjbL61iwg)++ および ++[各衛星の解析の準備](/2Q3bs6LQTlmidRxUzomJ6Q)++ が既に完了していることを前提とする。
すでにこのページを一読した人のためにScriptのみのページも用意した。
-> ++[NuSTARの解析(Scriptのみ)](/mzPYshyQSsCzrykukuAO4g)++

----

[TOC]

## はじめに
NuSTARはX線観測衛星の中で唯一無二の重要性があるのだが、それはhard X線(硬X線、目安として$10$keV以上)を高精度に観測できるからである。
NuSTARのホームページでも[$2$--$80$ keVの観測ができることがウリ](https://hackmd.io/HvYvSrQLQoqaRGDCMXz1dw?view#NuSTAR)であると謳われている。(実際には露出時間などの関係か$20$--$40$ keV程度以上は誤差が激しく、全部は解析に用いないことがある。)

**NuSTARには同じ性能のカメラが二つ(FPMA、FPMB)あるので、データ処理の最後に二つのカメラのデータを合算して統計の良いデータを得る。**

:::danger
このページでは観測データのディレクトリ構成として次のようなものを推奨し、前提としている。

work (観測データまとめ)
- NGC1313_X1 (天体名)
- data (観測データ)
- nustar (衛星名)
- 30XXXXXX (1つの観測データ)
- newton
- 06XXXXXX
:::

以下、具体的な操作手順を説明とともに記す。
説明のために分割したScriptのまとめを、各節ごとに「まとめ」として記載している。
基本的に「まとめ」の部分のコードブロックは、それぞれ`hea/work/<天体名>/data/<観測衛星>`でまとめてterminalにコピペできるようになっている。
##
### CALDBのダウンロード
++[各Mission(衛星)のCALDB一覧のページ](https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/caldb_supported_missions.html)++ からNuSTARの欄の++tar file++ をダウンロードし、CALDBディレクトリで展開する。

# 「Xspecの導入」でCALDBの設定は完了している前提。
# linkは最新のものに書き換える
cd $CALDB
wget http://heasarc.gsfc.nasa.gov/FTP/caldb/data/nustar/fpm/goodfiles_nustar_fpm.tar.gz -P $CALDB --no-check-certificate & pids="$!"
wait $pids && cd $CALDB && ls goodfiles* | xargs -n 1 -i tar xvf {}

:::info
nustar以外の衛星も上のリンク先にCALDBがあるものについては同様にしてCALDBが準備できる。
:::

:::info
`tar`による解凍では`*`(Wildcard)を使うことができない。そこで、上記では`ls goodfiles* | xargs -n 1 -i tar xvf {}`のように`xargs`を使っている。
簡単に説明すると`ls goodfiles*`による出力を(空白区切りで)1つずつ次のコマンドに(今の場合は`tar xvf`)に引き渡している。オプションの意味は以下の通り。
- `-n 1`: 1個ずつ引数に渡す。
- `-i`: xargsに渡された引数の置く場所を`{}`により明示的に表示することにする。
:::

:::warning
tarファイルをインストールしても何故かcaldb.indxだけ空のシンボリックリンクが作成される。(ターミナルで青い文字になる)ので、caldb.indxだけもう一度ダウンロードする。
https://heasarc.gsfc.nasa.gov/FTP/caldb/data/nustar/fpm/　

:::

## 1. nupipeline
まずは`nupipeline`で観測データに対して簡単な処理をする。
`work/<天体名>/data/nustar`でterminalに「まとめ」のscriptをコピペする。

### まとめ1

My_Nustar_D=$(pwd)
#cd $My_Newton_D
cd $My_Nustar_D
obs_dirs=($(find . -maxdepth 1 -type d -printf "%P\n" | grep ^[0-9]))
for dir in ${obs_dirs[@]}; do

    My_Nustar_ID=$dir
    My_Nustar_Dir=$My_Nustar_D/$My_Nustar_ID
    cd $My_Nustar_Dir

    rm -r -f $My_Nustar_Dir/out
    nupipeline indir=$My_Nustar_Dir \
        steminputs="nu$My_Nustar_ID" \
        outdir="$My_Nustar_Dir/out" saacalc="1" \
        saamode="optimized" tentacle="yes"
    #失敗したら`rm -r $My_Nustar_Dir/out`
done
cd $My_Nustar_D

:::warning
画像の一次処理らしい。実行すると10分くらいかかった。

:::

:::info
`My_Nustar_D=$(pwd)`で現在のpathを取得。
`obs_dirs=($(find . -maxdepth 1 -type d -printf "%P\n" | grep ^[0-9]))`で配列として観測データのディレクトリの一覧を取得。観測データのディレクトリ名は観測IDなので必ず数字から始まることで識別している。つまり、このnustarディレクトリに他のディレクトリを作成する場合は、数字から始まる名前は避ける。

bashのforループでディレクトリ名(=観測ID)を取得して`nupipeline`を実行している。
出力ディレクトリのoutが事前にあるとエラーになるので、直前に`rm`している。引数`-f`をつけることで、仮にoutがない場合もエラーを出さないようにしている。
:::

:::info
`nupipeline`の引数の`saacalc="1" saamode="optimized" tentacle="yes"`はフレアがあった場合に除去する用のオプション。

フレアチェック
: ++[NuSTARのフィルタチェック](http://www.srl.caltech.edu/NuSTAR_Public/NuSTAROperationSite/SAA_Filtering/SAA_Filter.php)++ のサイトで観測IDからフレア(明るすぎるBackground)が映っていないか確認できる。普通は自動フィルタできちんと除去されているし、`nupipeline`でも除去される。
:::

## 2. ds9で領域指定

こちらもnustarディレクトリでterminalにコピペする。NuSTARにはカメラが二つ(FPMA、FPMB)あるので一つの観測IDにつき二度ds9を起動する。また、一度のds9につき~~光源(source)とbackgroundの二つの領域ファイルを~~ **光源(source)とBackground(PropertyでBackgroundを指定)を指定する2つの円が含まれる1つの領域ファイルを保存する**。

### 具体的な領域指定の手順 (光源とBackgroundを含む1つの領域ファイルを保存する場合)
0. **事前に2つの円(点線の円と実線の円)を含む領域ファイルが読み込まれている場合**、それらを調整することで「1.-4.」を省略できる。
1. 画面中央のように光源天体が明らかに分かる場合はそこをクリックすると領域(黄緑色の円)とそのパラメータが載ったウィンドウが現れる。
2. 半径などに目的に沿う値を入力し、適用をクリックする。(先行研究や指導教官、先輩などに従えばよいが、とりあえず光源がカバーできていて、別の光源にかかっていなければよい。)
3. 続けて、Background領域のどこかをクリックして新たに領域(円)を指定する。(天体にかかっておらず、さして特徴のないあたりがよい。)
4. Background領域の円についても半径などを適当に調整する。そして、**追加ウィンドウの「属性(Property)」をクリックし、一番下の欄の「バックグラウンド領域」** をクリックする。すると、画像中の黄緑色の円が**実線から点線に切り替わる**。
5. その後、画面中央のタブの領域をクリックし、その下の段に現れる保存をクリックする。領域名は観測装置に合わせた名前にするとよい。(このページでは`fpmA.reg`、`fpmB.reg`のようにする。)
6. ds9のウィンドウを閉じると次のウィンドウが起動するので、同様の操作を延々続ける。**2つ目以降は前回の領域ファイルのコピーがあらかじめ読み込まれている**ので、それを調整して保存すればよい。

<img src="https://i.imgur.com/GDxOB7F.png" width="70%"></img>

<details><summary>具体的な領域指定の手順 (光源とBackgroundの2つの領域ファイルを個別に保存する場合)</summary>

### 具体的な領域指定の手順 (光源とBackgroundの2つの領域ファイルを個別に保存する場合)
1. 画面中央のように光源天体が明らかに分かる場合はそこをクリックすると領域(黄緑色の円)とそのパラメータが載ったウィンドウが現れる。
2. 半径などに目的に沿う値を入力し、適用をクリックする。(先行研究や指導教官、先輩などに従えばよいが、とりあえず光源がカバーできていて、別の光源にかかっていなければよい。)
3. その後、画面中央のタブの領域をクリックし、その下の段に現れる保存をクリックする。(このページでは)光源についてはディレクトリは変更せずに`sourceA.reg`で名前を付けて保存する。(FPMBの画像ならA->Bにする。)
4. さらに今ある円を(特にこだわりはなく)特徴のないBackground領域のどこかに移動させ、必要があれば半径を変更して`bkgA.reg`で保存する。(FPMBの画像ならA->Bにする。)
5. ds9のウィンドウを閉じると次のウィンドウが起動するので、同様の操作を延々続ける。**正直、領域指定が一番面倒。**

<img src="https://i.imgur.com/1meeEI1.png" width="70%"></img>

</details>

---

はじめて作業を行う場合など、forループを用いない処理の場合、はじめに観測IDなどをbash変数に保存する。

#My_Nustar_D=$(pwd)
cd $My_Nustar_D
obs_dirs=($(find . -maxdepth 1 -type d -printf "%P\n" | grep ^[0-9]))

dir_ind=0

My_Nustar_ID=${obs_dirs[$dir_ind]}
My_Nustar_Dir=$My_Nustar_D/$My_Nustar_ID

#cd $My_Nustar_DIr/fit
cd $My_Nustar_Dir/fit

### ds9起動

evt_file=nu${My_Nustar_ID}A01_cl.evt 
# evt_file=nu${My_Nustar_ID}B01_cl.evt

ds9 $evt_file -scale log \
-cmap bb -mode region

:::info
`ds9`の引数はそれぞれ次の意味を持つ。これらは(筆者が)ds9起動後によくつけていたオプションである。
- `-scale log` : logスケールで表示する。
- `-cmap bb` : カラーマップを黄～橙にする。(見やすい気がする)
- `-mode region` : 領域指定モードにする。
:::

変数の内容を更新して、別の観測IDに対しても同様に行う。

dir_ind=$(($dir_ind + 1))
if [[ ${#obs_dirs[@]} -le $dir_ind ]]; then
    My_Nustar_ID=${obs_dirs[$dir_ind]}
    My_Nustar_Dir=$My_Nustar_D/$My_Nustar_ID
    cd $My_Nustar_Dir/fit
else
    echo "The part is over."
fi

### まとめ2

#export My_Nustar_D=$(pwd)
cd $My_Nustar_D
obs_dirs=($(find . -maxdepth 1 -type d -printf "%P\n" | grep ^[0-9]))
for dir in ${obs_dirs[@]}; do

    My_Nustar_ID=$dir
    My_Nustar_Dir=$My_Nustar_D/$My_Nustar_ID
    cd $My_Nustar_Dir/out

    evt_file=nu${My_Nustar_ID}A01_cl.evt
    
    if [[ ! -f ${My_Nustar_D}/saved.reg ]]; then
        # saved.regが存在しないなら、新たに作成する
        declare -A tmp_dict=(["RA_OBJ"]="0" ["DEC_OBJ"]="0")
        for key in ${!tmp_dict[@]}; do
            # fits headerから座標を読み取る
            tmp_dict[$key]=$(fkeyprint infile="${evt_file}+1" keynam="${key}" |\
                grep "${key}\s*=" |\
                sed -r -n "s/^.*${key}\s*=\s*(.*)\s*\/.*$/\1/p")
        done
        ra=$(echo ${tmp_dict[RA_OBJ]} |\
            sed -r "s/E([\+\-]?[0-9]+)/*10^\1/" |\
            sed -r "s/10\^\+?(-?)0*([0-9]+)/10^(\1\2)/" | bc)
        dec=$(echo ${tmp_dict[DEC_OBJ]} |\
            sed -r "s/E([\+\-]?[0-9]+)/*10^\1/" |\
            sed -r "s/10\^\+?(-?)0*([0-9]+)/10^(\1\2)/" | bc)
        # background circleはとりあえず0.05 degずつずらした点
        ra_bkg=$(echo "$ra + 0.05 " | bc )
        dec_bkg=$(echo "$dec + 0.05 " | bc )
        # 半径はとりあえず0.026 deg = 100 arcsec
        ds9 $evt_file  \
            -regions system fk5 \
            -regions command "fk5; circle $ra $dec 0.026 # source" \
            -regions command "fk5; circle $ra_bkg $dec_bkg 0.026 # background" \
            -regions save $My_Nustar_D/saved.reg -exit
    fi
    cp ${My_Nustar_D}/saved.reg fpmA.reg -f

    ds9 nu${My_Nustar_ID}A01_cl.evt \
        -scale log -cmap bb -mode region \
        -regions load fpmA.reg
    ## adjust fpmA.reg

    cp fpmA.reg fpmB.reg -f
    ds9 nu${My_Nustar_ID}B01_cl.evt \
        -scale log -cmap bb -mode region \
        -region load fpmB.reg
    ## adjust fpmB.reg
    
    cp fpmB.reg ${My_Nustar_D}/saved.reg -f
    
    ## generate src[AB].reg and bkg[AB].reg from fpm[AB].reg
    for reg_file in $(find . -name "fpm[AB].reg"); do
        cam=$(echo ${reg_file} | sed -r -n "s/^.*fpm(A|B)\.reg$/\1/p")
        cat ${reg_file} | grep -v -E "^circle.*# background" > src${cam}.reg
        cat ${reg_file} | grep -v -E "^circle.*\)$" > bkg${cam}.reg
    done
done

cd $My_Nustar_D

:::warning
srcは半径4.5分角、bkgは半径1分角にした。suzakuの時と違って変にsaved.regとか作らなくてもうまくいった(と思う)。
今は/hea/work/mrk3/data/nustar/60002048004/out/にfpmA.regとfpmB.reg等がある。
→/hea/work/mrk3/data/nustar/60002048004/event_clに移動するべき？
:::

## 3. nuproductsでいろいろ生成
NuSTARの解析が楽な理由の最たるものは、この`nuproducts`でスペクトル解析に必要なファイルのほぼすべて(spectrum, lc, rmf, arf)を生成できるからである。**かなり時間がかかる。**

### まとめ3

#export My_Nustar_D=$(pwd)
cd $My_Nustar_D
obs_dirs=($(find . -maxdepth 1 -type d -printf "%P\n" | grep ^[0-9]))
for dir in ${obs_dirs[@]}; do

    My_Nustar_ID=$dir
    My_Nustar_Dir=$My_Nustar_D/$My_Nustar_ID

    cd $My_Nustar_Dir

    mkdir $My_Nustar_Dir/fit -p
    rm $My_Nustar_Dir/fit/* -f
    nuproducts \
        srcregionfile=$My_Nustar_Dir/out/srcA.reg \
        bkgregionfile=$My_Nustar_Dir/out/bkgA.reg \
        indir=$My_Nustar_Dir/out \
        outdir=$My_Nustar_Dir/fit \
        instrument=FPMA \
        steminputs=nu${My_Nustar_ID} \
        bkgextract=yes \
        clobber=yes

    nuproducts \
        srcregionfile=$My_Nustar_Dir/out/srcB.reg \
        bkgregionfile=$My_Nustar_Dir/out/bkgB.reg \
        indir=$My_Nustar_Dir/out \
        outdir=$My_Nustar_Dir/fit \
        instrument=FPMB \
        steminputs=nu${My_Nustar_ID} \
        bkgextract=yes \
        clobber=yes

done
cd $My_Nustar_D

:::warning

-------------------- nuproducts  error -----------------------
nuproducts_0.3.3: Unable to copy file: '356789_bk.pha' in '/home/kanta/hea/work/mrk3/data/nustar/60002048004/fit/nu60002048004B01_bk.pha'. No such file or directory
--------------------------------------------------------------

が出てきた。が、めげずにスクリプトを打ってはエラーを出すのを3回くらい繰り返していると急に成功した。
:::

:::info
light-curve

#export My_Nustar_D=$(pwd)
cd $My_Nustar_D
obs_dirs=($(find . -maxdepth 1 -type d -printf "%P\n" | grep ^[0-9]))
for dir in ${obs_dirs[@]}; do

    My_Nustar_ID=$dir
    My_Nustar_Dir=$My_Nustar_D/$My_Nustar_ID

    cd $My_Nustar_Dir

    mkdir $My_Nustar_Dir/fit -p
    rm $My_Nustar_Dir/fit/* -f
    nuproducts \
        srcregionfile=$My_Nustar_Dir/out/srcA.reg \
        bkgregionfile=$My_Nustar_Dir/out/bkgA.reg \
        indir=$My_Nustar_Dir/out \
        outdir=$My_Nustar_Dir/fit \
        instrument=FPMA \
        steminputs=nu${My_Nustar_ID} \
        binsize=1000 \
        pilow=800 \
        pihigh=4095\
        bkgextract=yes \
        clobber=yes

    nuproducts \
        srcregionfile=$My_Nustar_Dir/out/srcB.reg \
        bkgregionfile=$My_Nustar_Dir/out/bkgB.reg \
        indir=$My_Nustar_Dir/out \
        outdir=$My_Nustar_Dir/fit \
        instrument=FPMB \
        steminputs=nu${My_Nustar_ID} \
        binsize=1000 \
        pilow=800 \
        pihigh=4095 \
        bkgextract=yes \
        clobber=yes

done
cd $My_Nustar_D

:::

## 4. ひもづけとgrppha

はじめて作業を行う場合など、forループを用いない処理の場合、はじめに観測IDなどをbash変数に保存する。

#My_Nustar_D=$(pwd)
cd $My_Nustar_D
obs_dirs=($(find . -maxdepth 1 -type d -printf "%P\n" | grep ^[0-9]))

dir_ind=0

My_Nustar_ID=${obs_dirs[$dir_ind]}
My_Nustar_Dir=$My_Nustar_D/$My_Nustar_ID

cd $My_Nustar_Dir/fit

### addascaspec
同じ観測の光源データ、Backgroundデータ、およびarf、rmfファイルを紐づけする操作。(arf、rmfの説明はSuzakuの解析あたりででも説明予定)

cat <<EOF > tmp_fi.add
nu${My_Nustar_ID}A01_sr.pha nu${My_Nustar_ID}B01_sr.pha
nu${My_Nustar_ID}A01_bk.pha nu${My_Nustar_ID}B01_bk.pha
nu${My_Nustar_ID}A01_sr.arf nu${My_Nustar_ID}B01_sr.arf
nu${My_Nustar_ID}A01_sr.rmf nu${My_Nustar_ID}B01_sr.rmf
EOF

rm AB_${My_Nustar_ID}_nongrp.fits \
    AB_${My_Nustar_ID}_rsp.fits \
    AB_${My_Nustar_ID}_bkg.fits -f
addascaspec tmp_fi.add \
    AB_${My_Nustar_ID}_nongrp.fits \
    AB_${My_Nustar_ID}_rsp.fits \
    AB_${My_Nustar_ID}_bkg.fits

:::warning
実行結果は以下の通り。

The fourth parameter AB_60002048004_bkg.fits is interpreted as outbgd.
Number of spectral files to add = 2
### Set no 1 ###
Spec:nu60002048004A01_sr.pha
  Exposure time:  3.339e+04
  Backscal     :  3.710e-02
Bgd :nu60002048004A01_bk.pha
  Exposure time:  3.339e+04
  Backscal     :  1.874e-03
ARF :nu60002048004A01_sr.arf
RMF :nu60002048004A01_sr.rmf

### Set no 2 ###
Spec:nu60002048004B01_sr.pha
  Exposure time:  3.341e+04
  Backscal     :  3.682e-02
Bgd :nu60002048004B01_bk.pha
  Exposure time:  3.341e+04
  Backscal     :  1.875e-03
ARF :nu60002048004B01_sr.arf
RMF :nu60002048004B01_sr.rmf

mathpha expr="@/tmp/14721.mathpha_expr" units=C outfil="AB_60002048004_nongrp.fits" exposure=CALC \
properr=no areascal=NULL auxfiles=NONE backfile=NONE arfile=NONE rmfile=NONE \
backscal=1.0E0 ncomments=1 comment1="Produced by addascaspec version 1.30"
 ......... processing file:
 ......... nu60002048004A01_sr.pha
 ......... processing file:
 ......... nu60002048004B01_sr.pha
 WARNING - populate 5.3.0: Instrument mismatch
 ...... data from different mission/inst/filter than prev datasets
 WARNING - populate 5.3.0: Instrumental discrepancies detected
 ...... proceeding, but maybe you should check your i/p files
 ... Following keywrds set to null values in o/p file:
 ...... TELESCOP, INSTRUME, FILTER
 ...  ... performing algebra in units of COUNTS
 ... wrote the PHA data Extension
 ** mathpha 4.7.0 completed successfully
fparkey value=1.0 fitsfile=AB_60002048004_nongrp.fits keyword=BACKSCAL add=yes
fparkey value="AB_60002048004_bkg.fits" fitsfile=AB_60002048004_nongrp.fits keyword=BACKFILE add=yes

mathpha expr="@/tmp/14721.mathpha_expr" units=C outfil="AB_60002048004_bkg.fits" exposure=CALC \
properr=no areascal=NULL auxfiles=NONE backfile=NONE arfile=NONE  rmfile=NONE \
backscal=5.070e-02 ncomments=1 comment1="Produced by addascaspec version 1.30"
 ......... processing file:
 ......... nu60002048004A01_bk.pha
 ......... processing file:
 ......... nu60002048004B01_bk.pha
 WARNING - populate 5.3.0: Instrument mismatch
 ...... data from different mission/inst/filter than prev datasets
 WARNING - populate 5.3.0: Instrumental discrepancies detected
 ...... proceeding, but maybe you should check your i/p files
 ... Following keywrds set to null values in o/p file:
 ...... TELESCOP, INSTRUME, FILTER
 ...  ... performing algebra in units of COUNTS
 ... wrote the PHA data Extension
 ** mathpha 4.7.0 completed successfully
fparkey value=5.070e-02 fitsfile=AB_60002048004_bkg.fits+1 keyword=BACKSCAL add=yes
fparkey value="AB_60002048004_rsp.fits" fitsfile=AB_60002048004_nongrp.fits keyword=RESPFILE add=yes
marfrmf rmfil=nu60002048004A01_sr.rmf arfil=nu60002048004A01_sr.arf outfil=/tmp/147210.rsp qoverride=no
 ** marfrmf 3.3.1
 ...  HDUCLAS3 in o/p RSP_MATRIX set to FULL
keyword not found in header
 RDARF1 1.1.0   WARNING:  reading FILTER
keyword not found in header
 RDARF1 1.1.0   WARNING:  reading DETNAM
 WARNING - do_marfrmf 4.4.0:  Following conflicts found:
 ... DETNAM   from RMF  = DET0
 ... DETNAM   from ARF  = NONE
 ...... values from RMF will be used in o/p
 ** marfrmf 3.3.1 completed successfully
marfrmf rmfil=nu60002048004B01_sr.rmf arfil=nu60002048004B01_sr.arf outfil=/tmp/147211.rsp qoverride=no
 ** marfrmf 3.3.1
 ...  HDUCLAS3 in o/p RSP_MATRIX set to FULL
keyword not found in header
 RDARF1 1.1.0   WARNING:  reading FILTER
keyword not found in header
 RDARF1 1.1.0   WARNING:  reading DETNAM
 WARNING - do_marfrmf 4.4.0:  Following conflicts found:
 ... DETNAM   from RMF  = DET0
 ... DETNAM   from ARF  = NONE
 ...... values from RMF will be used in o/p
 ** marfrmf 3.3.1 completed successfully

addrmf @/tmp/14721.addrmf_list rmffile=AB_60002048004_rsp.fits
addrmf v1.23
Summing ...
 5.000E-01 * /tmp/147210.rsp
 5.000E-01 * /tmp/147211.rsp
 RMF #    1 : /tmp/147210.rsp            0.50000
              NuSTAR     FPMA       DET0       NONE       PI
 RMF #    2 : /tmp/147211.rsp            0.50000
              NuSTAR     FPMB       DET0       NONE       PI
 ... using rdrmf5 1.3.3
 ... read all the keywords
 ... reading the data
 ... successfully read RSP_MATRIX data
 ... using rdrmf5 1.3.3
 ... read all the keywords
 ... reading the data
 ... successfully read RSP_MATRIX data
 ...  using wtrmf5 1.5.0
 ... new extension created
 ...  Max # elements in F_CHAN & N_CHAN arrays =        7
 ...  Max # elements in MATRIX array =     2006
 ... Using Variable length array for MATRIX column
 ... written the extension header keywords
 ... written the OGIP required keywords
 ... written the history keywords
 ... written the comment keywords
 ... defined the data structure
 ... successfully written RSP_MATRIX data

rm -f /tmp/14721.addrmf_list
rm -f /tmp/147210.rsp
rm -f /tmp/147211.rsp
rm -f /tmp/14721.mathpha_expr

:::
### grppha
これはbinまとめと呼ばれる作業であり、観測データを扱う際には必須な操作だ。扱うエネルギー帯や観測天体にもよるがbinの数に対して入力数が十分に多くない(**まばらな**)場合、統計的に扱えなくなってしまう。
そこで各binをいくつかずつ束ねて一つのbinと見なす(++binまとめ++)ことで、統計的なふるまいを扱うことができるようになる。
今回のScriptではカメラAおよびBに対しては入力が50以上になるように、A+B(AとBの合算)では100以上になるようにbinをまとめさせている。このまとめの単位は解析ごとに適宜判断する。

gnum=50
grp_name=AB_${My_Nustar_ID}_grp${gnum}.fits
rm ${grp_name} -f
grppha infile=AB_${My_Nustar_ID}_nongrp.fits outfile=${grp_name}
group min 50
exit

変数の内容を更新して、別の観測IDに対しても同様に行う。

dir_ind=$(($dir_ind + 1))
if [[ ${#obs_dirs[@]} -le $dir_ind ]]; then
    My_Nustar_ID=${obs_dirs[$dir_ind]}
    My_Nustar_Dir=$My_Nustar_D/$My_Nustar_ID
    cd $My_Nustar_Dir/fit
else
    echo "The part is over."
fi

### まとめ4

## ひもづけとgrppha
gnum=50
echo ${My_Nustar_D:=$(pwd)} # 未定義時に代入
cd $My_Nustar_D
obs_dirs=($(find . -maxdepth 1 -type d -printf "%P\n" | grep ^[0-9]))
for dir in ${obs_dirs[@]}; do

My_Nustar_ID=$dir
My_Nustar_Dir=$My_Nustar_D/$My_Nustar_ID

cd $My_Nustar_Dir/fit/

cat <<EOF > tmp_fi.add
nu${My_Nustar_ID}A01_sr.pha nu${My_Nustar_ID}B01_sr.pha
nu${My_Nustar_ID}A01_bk.pha nu${My_Nustar_ID}B01_bk.pha
nu${My_Nustar_ID}A01_sr.arf nu${My_Nustar_ID}B01_sr.arf
nu${My_Nustar_ID}A01_sr.rmf nu${My_Nustar_ID}B01_sr.rmf
EOF

rm AB_${My_Nustar_ID}_nongrp.fits \
    AB_${My_Nustar_ID}_rsp.fits \
    AB_${My_Nustar_ID}_bkg.fits -f
addascaspec tmp_fi.add \
    AB_${My_Nustar_ID}_nongrp.fits \
    AB_${My_Nustar_ID}_rsp.fits \
    AB_${My_Nustar_ID}_bkg.fits

grp_name=AB_${My_Nustar_ID}_grp${gnum}.fits

## edit header for spectrum file
oldName=nu${My_Nustar_ID}A01_sr.pha
newName=$nongrp_name

### same values
cp_keys=(TELESCOP OBS_ID TARG_ID OBJECT RA_OBJ \
DEC_OBJ RA_NOM DEC_NOM RA_PNT DEC_PNT PA_PNT \
EQUINOX RADECSYS TASSIGN TIMESYS MJDREFI MJDREFF \
TIMEREF CLOCKAPP TIMEUNIT TSTOP DATE-OBS DATE-END \
ORIGIN DATETLM TLM2FITS SOFTVER CALDBVER USER \
TIMEZERO DEADAPP TSORTKEY NUFLAG ABERRAT FOLOWSUN \
DSTYP1 DSVAL1 NUPSDOUT DEPTHCUT OBSMODE HDUNAME \
AXLEN1 AXLEN2 CONTENT WMREBIN OBS-MODE SKYBIN \
PIXSIZE WMAPFIX DSTYP2 DSREF2 DSVAL2 CTYPE1 DRPIX1 \
CRVAL1 CDELT1 DDELT1 CTYPE2 DRPIX2 CRVAL2 CDELT2 \
DDELT2 WCSNAMEP WCSTY1P LTM1_1 CTYPE1P CRPIX1P \
CDELT1P WCSTY2P LTM2_2 CTYPE2P CRPIX2P CDELT2P \
OPTIC1 OPTIC2 HBBOX1 HBBOX2 REFXCTYP REFXCRPX \
REFXCRVL REFXCDLT REFYCTYP REFYCRPX REFYCRVL REFYCDLT)

### near values
cp_keys2=(INSTRUME TSTART TELAPSE ONTIME LIVETIME \
MJD-OBS FILIN001 DEADC NPIXSOU CRPIX1 CRPIX2 LTV1 \
CRVAL1P LTV2 CRVAL2P BBOX1 BBOX2 X-OFFSET \
Y-OFFSET TOTCTS)

declare -A tr_keys=(
    ["BACKFILE"]="AB_${My_Nustar_ID}_bkg.fits"
    ["RESPFILE"]="AB_${My_Nustar_ID}_rsp.fits"
)

for key in ${cp_keys[@]} ${cp_keys2[@]}; do
    orig_val=$(fkeyprint infile="${oldName}+0" keynam="${key}" |\
        grep "${key}\s*=" |\
        sed -r -n "s/^.*${key}\s*=\s*(.*)\s*\/.*$/\1/p")

    tr_keys[$key]="${orig_val}"
done

for key in ${!tr_keys[@]}; do
    fparkey value="${tr_keys[$key]}" \
        fitsfile=${newName}+1 \
        keyword="${key}" add=yes
done

## edit header for bkg file
oldName=nu${My_Nustar_ID}A01_bk.pha
newName=AB_${My_Nustar_ID}_bkg.fits

### same values
cp_keys=(TELESCOP OBS_ID TARG_ID OBJECT RA_OBJ \
DEC_OBJ RA_NOM DEC_NOM RA_PNT DEC_PNT PA_PNT \
EQUINOX RADECSYS TASSIGN TIMESYS MJDREFI MJDREFF \
TIMEREF CLOCKAPP TIMEUNIT TSTART TSTOP TELAPSE \
DATE-OBS DATE-END ORIGIN CREATOR DATETLM TLM2FITS \
SOFTVER CALDBVER MJD-OBS USER FILIN001 TIMEZERO \
DEADAPP TSORTKEY NUFLAG ABERRAT FOLOWSUN DSTYP1 \
DSVAL1 NPIXSOU NUPSDOUT DEPTHCUT OBSMODE HDUNAME \
AXLEN1 AXLEN2 CONTENT WMREBIN OBS-MODE SKYBIN \
PIXSIZE WMAPFIX DSTYP2 DSREF2 DSVAL2 CTYPE1 CRPIX1 \
DRPIX1 CRVAL1 CDELT1 DDELT1 CTYPE2 CRPIX2 DRPIX2 \
CRVAL2 CDELT2 DDELT2 WCSNAMEP WCSTY1P LTV1 LTM1_1 \
CTYPE1P CRPIX1P CRVAL1P CDELT1P WCSTY2P LTV2 \
LTM2_2 CTYPE2P CRPIX2P CRVAL2P CDELT2P OPTIC1 \
OPTIC2 BBOX1 BBOX2 HBBOX1 HBBOX2 X-OFFSET Y-OFFSET \
REFXCTYP REFXCRPX REFXCRVL REFXCDLT REFYCTYP \
REFYCRPX REFYCRVL REFYCDLT)

### near values
cp_keys2=(INSTRUME DATE ONTIME LIVETIME DEADC)

declare -A tr_keys=()

for key in ${cp_keys[@]} ${cp_keys2[@]}; do
    orig_val=$(fkeyprint infile="${oldName}+0" keynam="${key}" |\
        grep "${key}\s*=" |\
        sed -r -n "s/^.*${key}\s*=\s*(.*)\s*\/.*$/\1/p")

    tr_keys[$key]="${orig_val}"
done

for key in ${!tr_keys[@]}; do
    fparkey value="${tr_keys[$key]}" \
        fitsfile=${newName}+1 \
        keyword="${key}" add=yes

rm ${grp_name} -f
cat << EOF | bash
grppha infile=AB_${My_Nustar_ID}_nongrp.fits outfile=${grp_name} clobber=true
group min ${gnum}
exit !${grp_name}
EOF

done
done
cd $My_Nustar_D

:::info
ここの「まとめ」ではnongrpファイルのheaderも編集している。
これは、改めて後にbinまとめをする場合に対応している。

また、addascaspecで相当のheader情報が抜けてしまっているので、元ファイルのheaderから、問題のない情報のみコピーしている。
:::
:::warning
bin50だと今回は足りなかったようでxspecの画像が下のようになった。
bin500だと以下。
下振れがなくなってきた。
:::

## 5. fitディレクトリにまとめ

解析に必要なファイルのみを別のディレクトリにまとめる。特に複数の衛星のデータをまとめて解析することが多いので、nustarディレクトリと同じ階層にもfitディレクトリを作って(`$My_Nustar_D/../fit/`)そこにシンボリックリンクの形で各衛星の解析用のファイルを集める。

移行の作業(=スペクトル解析)は`$My_Nustar_D/../fit/`で行う。

### まとめ5

## fitsディレクトリにまとめ
echo ${My_Nustar_D:=$(pwd)} 
cd $My_Nustar_D
obs_dirs=($(find . -maxdepth 1 -type d -printf "%P\n" | grep ^[0-9]))
mkdir -p $My_Nustar_D/fit $My_Nustar_D/../fit/
for dir in ${obs_dirs[@]};do
    cp $dir/fit/AB_* fit/
done
find $My_Nustar_D/fit/ -name "AB_*.*" -type f  -printf "%f\n" |\
    xargs -n 1 -i rm -f $My_Nustar_D/../fit/{}
ln -s $My_Nustar_D/fit/AB_*.* $My_Nustar_D/../fit/
cd $My_Nustar_D

:::warning
binまとめ成功？のbin500のデータだけシンボリックリンクができていない。
mrk3/data/fitのディレクトリで

ln -s ../nustar/fit/AB_60002048004_grp500.fits

を実行するとできた。
:::