MadGraph5_aMC@NLO, MadAnalysis5

(under construction)

MadGraphでscanコマンドを用いて、例えば粒子の質量を変化させると、その都度対応する断面積のデータがテキストファイルに書き出される。このデータをgnuplot等の描画ソフトでグラフに描くことができるが、ここではpythonを使ってグラフを作成する方法を紹介する。

次のような2種類の(x,y)データを用意する。それぞれ模型内の粒子質量に対する断面積の値をまとめたもので、二種類のデータは模型のあるパラメータの値を変えた場合に対応している：

file名：scan_mx10[-_scan_10].txt

#run_name            mass#5000001         cross(pb)

mx10                 2.000000e+02         1.149446e+03

mx10_scan_02         4.000000e+02         1.365383e+02

mx10_scan_03         6.000000e+02         3.545218e+01

mx10_scan_04         8.000000e+02         1.278913e+01

mx10_scan_05         1.000000e+03         5.467218e+00

mx10_scan_06         1.200000e+03         2.653384e+00

mx10_scan_07         1.400000e+03         1.391719e+00

mx10_scan_08         1.600000e+03         7.551432e-01

mx10_scan_09         1.800000e+03         4.382577e-01

mx10_scan_10         2.000000e+03         2.609864e-01

file名：less scan_mx_05my1m10\[-_scan_10\].txt

#run_name            mass#5000001         mass#5000521         cross(pb)

mx_05my1m10          2.000000e+02         9.000000e+01         8.261480e+02

mx_05my1m10_scan_02  4.000000e+02         1.900000e+02         7.624530e+01

mx_05my1m10_scan_03  6.000000e+02         2.900000e+02         1.640644e+01

mx_05my1m10_scan_04  8.000000e+02         3.900000e+02         5.252261e+00

mx_05my1m10_scan_05  1.000000e+03         4.900000e+02         2.031126e+00

mx_05my1m10_scan_06  1.200000e+03         5.900000e+02         9.142628e-01

mx_05my1m10_scan_07  1.400000e+03         6.900000e+02         4.399481e-01

mx_05my1m10_scan_08  1.600000e+03         7.900000e+02         2.262258e-01

mx_05my1m10_scan_09  1.800000e+03         8.900000e+02         1.231878e-01

mx_05my1m10_scan_10  2.000000e+03         9.900000e+02         6.894974e-02

この2つのデータファイルと同じディレクトリに例えば"plot.py"というファイル名で次のようなコードを書く（馬渡氏提供）：

# matplotlib example file by K. Mawatari (2017.11.28)

# type 'python plot.py'

import numpy as np

import matplotlib.pyplot as plt

# import data files

data1 = np.loadtxt('scan_mx10[-_scan_10].txt', skiprows=1, usecols=(1,2)) # 1D plot

data2 = np.loadtxt('scan_mx_05my1m10[-_scan_10].txt', skiprows=1, usecols=(1,2,3)) # 2D plot

# define the columns

my1=data1[:,0]

xsec1=data1[:,1]

xsec2=data2[:,2]

# define the variables for the x and y axes

plt.plot(my1,xsec1,'-o',label='$m_X=10$')

plt.plot(my1,xsec2,'--o',label='$m_X=m_Y/2-10$')

# set up the plots

xlab=r'$m_{Y}$ [GeV]'

ylab=r'$\sigma(pp\to XXj)$ [pb]'

xlimrange=[0,2000]

ylimrange=[0.01,10000]

savepdf='test.pdf'

plt.yscale('log')

plt.xlim(xlimrange)

plt.ylim(ylimrange)

plt.title("Mono-j for simplified DM (s-channel spin-1)")

plt.xlabel(xlab, fontsize=18)

plt.ylabel(ylab, fontsize=18)

plt.legend()

plt.grid()

### make the plot

plt.savefig(savepdf, format='pdf', bbox_inches='tight')

plt.show()

ターミナルで

$ python plot.py

とすると次のようなグラフが表示される。グラフの見栄えを変えたりするのはググってみること。

scanコマンドを使って次のようなデータを作成する（ダウンロード：scan_maddm_[1-209].txt)

file名: scan_maddm_[1-209].txt

#run_name            mass#5000001         mass#5000521         omegah2              x_freezeout          sigmav_xf

1                    5.000000e+02         5.000000e+01         3.626097e-02         2.100000e+01         3.196735e-09

1                    5.000000e+02         1.000000e+02         7.296928e-03         2.400000e+01         1.704846e-08

2                    5.000000e+02         1.500000e+02         1.939136e-03         2.500000e+01         6.948923e-08

3                    5.000000e+02         2.000000e+02         3.031988e-04         2.800000e+01         5.295001e-07

4                    5.000000e+02         2.500000e+02         2.324942e-06         3.300000e+01         4.366346e-05

5                    5.000000e+02         3.000000e+02         2.479971e-04         2.800000e+01         4.773687e-07

6                    5.000000e+02         3.500000e+02         7.874220e-04         2.700000e+01         1.508970e-07

7                    5.000000e+02         4.000000e+02         1.524455e-03         2.600000e+01         7.709789e-08

8                    5.000000e+02         4.500000e+02         2.417971e-03         2.600000e+01         4.840366e-08

...

これを等高線表示させるコードは例えば次の通り（馬渡氏提供）：

# matplotlib example file for a contour plot by K. Mawatari (2017.11.29)

# type 'python plot_contour.py'

import numpy as np

import matplotlib.pyplot as plt

# import data files

data1 = np.loadtxt('scan_maddm_[1-209].txt', skiprows=1, usecols=(1,2,3)) # 3D data

# define the columns

#x=data1[:,0]

#y=data1[:,1]

x=np.linspace(500, 2500, 21)

y=np.linspace(50, 500, 10)

X,Y=np.meshgrid(x, y)

z=data1[:,2].reshape(21,10)

Z=z.transpose()

cont=plt.contour(X,Y,Z, levels=[0.001,0.01, 0.12, 1])

cont.clabel(fmt='%1.3f', fontsize=14)

# set up the plots

xlab=r'$m_Y$ [GeV]'

ylab=r'$m_X$ [GeV]'

xlimrange=[0,2500]

ylimrange=[0,500]

savepdf='test_contour.pdf'

#plt.yscale('log')

plt.xlim(xlimrange)

plt.ylim(ylimrange)

plt.title("Relic density for simplified DM (s-channel spin-1)")

plt.xlabel(xlab, fontsize=18)

plt.ylabel(ylab, fontsize=18)

#plt.legend()

plt.grid()

### make the plot

plt.savefig(savepdf, format='pdf', bbox_inches='tight')

plt.show()

これを例えば"plot_contour.py"という名前のファイルに保存し、

$ python plot_contour.py

と打てば次のグラフが表示される：

Report abuse