ダブルΛハイパー核 NagaraイベントをROOTとPythonを使ってMCで発生させてみる

3体崩壊で、粒子が完全にランダムに崩壊したとすると、粒子1と粒子2の不変質量の2乗と、粒子1と粒子3の不変質量の2乗分布が一様になるように崩壊する。

ROOT のライブラリを使って、これを実際に発生してみる。B_ΞとΔB_ΛΛは0としている。

#include <iostream>
#include <fstream>
#include <TLorentzVector.h>
#include <TGenPhaseSpace.h>
void PhaseSpace()
{
   TLorentzVector target(0.0, 0.0, 0.0, 11174.864 * 0.001);
   TLorentzVector beam(0.0, 0.0, 0.0, 1321.71 * 0.001);
   TLorentzVector W = beam + target;
   //(Momentum, Energy units are Gev/C, GeV)
   Double_t masses[3] = {(3727.379 + (1115.683 - 3.12) * 2) * 0.001, 3727.379 * 0.001, 2808.921 * 0.001};
   TGenPhaseSpace event;
   event.SetDecay(W, 3, masses);

   std::ofstream ofs("nagara.txt");
   for (Int_t n = 0; n < 100000; n++)
   {
      Double_t weight = event.Generate();
      Double_t theta = 0;
      Double_t phi = 0;
      for (Int_t p = 0; p < 3; p++)
      {
         auto pParticle = event.GetDecay(p);
         //4元運動量を出力
         ofs << pParticle->E() * 1000 << " ";
         ofs << pParticle->Px() * 1000 << " ";
         ofs << pParticle->Py() * 1000 << " ";
         ofs << pParticle->Pz() * 1000 << " ";
      }
      auto p1 = *event.GetDecay(0) + *event.GetDecay(1);
      auto p2 = *event.GetDecay(0) + *event.GetDecay(2);
      //He6LLとHe4の不変質量の2乗
      ofs << p1.M2() << " ";
      //He6LLとtの不変質量の2乗
      ofs << p2.M2() << " ";
      //イベントの重み
      ofs << weight << " ";
      ofs << "\n";
   }
   ofs.close();
}

低速ではあるが、インタプリタで実行した。

root PhaseSpace.C -q -l

それぞれの不変質量の2乗をX軸、Y軸に取ったときの2Dヒストグラム。いわゆるDalitz plot

Weightあり

Weightなし

Pythonだと、Albert Puig Navarro氏が開発しているphasespaceモジュールが便利。

pip install git+https://github.com/zfit/phasespace

で最新版がインストールできる。テストコードは以下の通り。

# %%
import matplotlib.pyplot as plt
import phasespace
import tensorflow as tf
import numpy as np
# %%
with tf.Graph().as_default():
    tf.compat.v1.set_random_seed(1)

# %%
B0_MASS = 11174.864+1321.71
P1_MASS = 3727.379 + (1115.683 - 3.12) * 2
P2_MASS = 3727.379
P3_MASS = 2808.921

weights, particles = phasespace.nbody_decay(B0_MASS, [P1_MASS, P2_MASS, P3_MASS]).generate(n_events=100000)

# %%
vm01 = []
vm02 = []
for i, (p0, p1, p2) in enumerate(zip(particles["p_0"], particles["p_1"], particles["p_2"])):
    p01 = np.array(p0)+np.array(p1)
    p02 = np.array(p0)+np.array(p2)
    p01 = p01*p01  # 自身の内積
    p02 = p02*p02  # 自身の内積
    m01 = p01[3]-(p01[0]+p01[1]+p01[2])  # 粒子1と粒子2の不変質量
    m02 = p02[3]-(p02[0]+p02[1]+p02[2])  # 粒子1と粒子3の不変質量
    vm01.append(m01*1e-6)
    vm02.append(m02*1e-6)

plt.hist2d(vm01, vm02, weights=weights, bins=[20, 20])
plt.show()
plt.hist2d(vm01, vm02, bins=[20, 20])
plt.show()

毎回同じ結果を生成させるために、下記のコードでSeedを固定しているが、固定の必要がないなら除けば良い。

with tf.Graph().as_default():
    tf.compat.v1.set_random_seed(1)

参考にしたサンプルコード

root.cern

qiita.com