---
lang:ja-jp
###### tags:
---

# PyXspecでmatplotlibを利用する

[X線スペクトル解析](/wOsMOU-3QCC6rkHt3kdfuA)

## はじめに

2020年2月にHEASoft Ver.6.27がリリースされ、その中でPyXspecに新たに`Plot.addComp()`という関数が追加された。これはモデルの各要素(Xspecの`plot eeu`における点線)の座標に対応する。逆に今までは、このモデルの各要素の座標をpython側で取り出すことができなかった。


### module
以下のモジュールをimportする。
`gridspec`はmatplotlibのsubplotにおけるグラフの縦幅を調整するのに使用する。

from xspec import *
import matplotlib.pyplot import plt
from matplotlib import gridspec

## plotデータを取得する

### x座標, y座標, 誤差
これらは今までのPyXspecでも取得できた。観測データ点の座標および誤差である。

xs = Plot.x(plotGroup, plotWindow)
ys = Plot.y(plotGroup, plotWindow)
xe = Plot.xErr(plotGroup, plotWindow)
ye = Plot.yErr(plotGroup, plotWindow)

引数の`plotGroup`はXspecにおけるData Groupの番号に対応する。PyXspecではDataManagerやModelのindexと対応する。DataManagerやModelはそれぞれmanualでは`s[0], s[1]`や`m[0], m[1]`に収納される。

`plotWindow`は`plot eeu ratio`などで複数のグラフがプロットされた場合に、グラフを指定するパラメーターだ。qdpにおける`window 1`などの番号とも対応している。上から順に1, 2, ...と数える。

### モデル(総和)
これも今までのPyXspecでも取得できた。各モデル(component)の総和に対応する。

ys_model = Plot.model(plotGroup, plotWindow)

### component(モデル各要素)
今回のPyXspecで初めて取得できるようになった。自分で設定するモデルに応じて要素数は変化するが、ここではwhileのループと`try, except`で適当にすべて取得している。

ys_comps = []
comp_N = 1
while(True):
    try:
        ys_tmp = Plot.addComp(comp_N, plotGroup, plotWindow)
        comp_N += 1
        # execlude components with only 0
        if sum([1 for s in ys_tmp if s == 0]) == len(ys_tmp):
            continue
        ys_comps.append(ys_tmp)
    except:
        break

### まとめ
各座標データが存在しているかどうかは何を`plot`するかによる。ここではいくつかを例に挙げる。
なお同時に複数のグラフをプロットした場合、x座標は共通である。

|plot|x座標, 誤差 |y座標, 誤差|モデル(総和)|component|
|-|-|-|-|-|
|eeu|O|O|O|O|
|eem|O|X|O|O|
|ratio|O|O|X|X|
|del|O|O|X|X|

以上をまとめると、次のような関数になる。

def obtain_xys(plots=["eeu"]):
    plot_command = " ".join(plots)
    Plot(plot_command)
    xys_s = {}

    for plotWindow_tmp, plot_type in enumerate(plots):
        plotWindow = plotWindow_tmp+1
        xys = {}
        xs, ys, xe, ye, ys_model, ys_comps = [[]]*6
        for plotGroup in range(1, AllData.nGroups+1):
            xs = Plot.x(plotGroup, plotWindow)
            if plot_type in {"eeu", "ratio", "del"}:
                ys = Plot.y(plotGroup, plotWindow)
                xe = Plot.xErr(plotGroup, plotWindow)
                ye = Plot.yErr(plotGroup, plotWindow)

            if plot_type in {"eeu", "eem"}:
                ys_model = Plot.model(plotGroup, plotWindow)

            # obtain comps in models
            ys_comps = []
            comp_N = 1
            while(True):
                try:
                    ys_tmp = Plot.addComp(comp_N, plotGroup, plotWindow)
                    comp_N += 1
                    # execlude components with only 0
                    if sum([1 for s in ys_tmp if s == 0]) == len(ys_tmp):
                        continue
                    ys_comps.append(ys_tmp)
                except:
                    break
            xys[plotGroup] = {"xs": xs, "ys": ys, "xe": xe,
                              "ye": ye, "ys_model": ys_model, "ys_comps": ys_comps}
        xys_s[plot_type]=xys

    return xys_s

## グラフを作成する
上で定義した`obtain_xys()`を利用して、matplotlibでXspec風にグラフを作成する。
ここまでならmatplotlibを利用した利点はないが、グラフの調整をmatplotlibからできるというのは、考えるだけでも嬉しくなってくる人が多いことだろう。


def plot_xys(plots=["eeu"], xlim=[], ylims={}, colorlist=["r", "g", "b", "c", "m", "y", "k", "w"]):
    plots_diff = set(plots)-{"eeu", "eem", "ratio", "del"}
    if len(plots_diff) >= 1:
        print(f"{', '.join(list(plots_diff))} are not appropriate")

    fig = plt.figure()
    fig.patch.set_facecolor('white')

    xys_s = obtain_xys(plots=plots)
    subplots = []

    # set height ratios for sublots
    gs = gridspec.GridSpec(len(plots), 1, height_ratios=[
                           2 if s in ["eeu", "eem"] else 1 for s in plots])

    for gs_tmp, plot_type in zip(gs, plots):
        xys = xys_s[plot_type]
        # the fisrt subplot
        ax = plt.subplot(gs_tmp)
        subplots.append(ax)

        # set label
        plt.subplots_adjust(hspace=.0)
        if not gs_tmp == gs[-1]:
            plt.setp(ax.get_xticklabels(), visible=False)
        else:
            plt.xlabel("Energy (keV)")
        ylabel_dic = {"eeu": r"photons cm$^\mathdefault{-2}$ s$^\mathdefault{-1}$ keV$^\mathdefault{-1}$",
                      "eem": "$EF_\mathdefault{E}$ (keV$^2$)", "ratio": "ratio", "del": "residual"}
        plt.ylabel(ylabel_dic[plot_type])

        fig.align_labels()

        # set range
        plt.xscale("log")

        if not plot_type in ["ratio", "del"]:
            plt.yscale("log")

        # plot eeu
        for plotGroup, xys_tmp in xys.items():
            xs = xys_tmp["xs"]
            if plot_type in {"eeu", "ratio", "del"}:
                ys = xys_tmp["ys"]
                xe = xys_tmp["xe"]
                ye = xys_tmp["ye"]

            # plt.scatter(xs,ys)
            plt.errorbar(xs, ys, yerr=ye, xerr=xe, capsize=0, fmt="o", markersize=5,
                         ecolor=colorlist[plotGroup], markeredgecolor=colorlist[plotGroup], color="none")
            if plot_type in {"eeu", "eem"}:
                ys_model = xys_tmp["ys_model"]
                plt.plot(xs, ys_model, color=colorlist[plotGroup])

            ys_comps = xys_tmp["ys_comps"]
            for ys_comp in ys_comps:
                plt.plot(xs, ys_comp, linestyle="dotted",
                         color=colorlist[plotGroup])
            # plt.plot()
    return fig, subplots, xys_s