スーパーカミオカンデのイベントディスプレイを作ってみよう

東京大学宇宙線研究所の神岡宇宙素粒子研究施設より、スーパーカミオカンデ実験のシミュレーションデータが公開されました。Python の描画スクリプト付きなので、遊んでみたいと思います。 

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">スーパーカミオカンデのイベントディスプレイ作成用データの公開を始めました。ニュートリノや陽子崩壊のシミュレーションデータを元に、あなた独自のイベントディスプレイを作ってみませんか？データには光電子増倍管の位置座標と光の量や時間の情報が含まれています。<a href="https://t.co/YTqHhOTHAP">https://t.co/YTqHhOTHAP</a> <a href="https://t.co/PZVXqA0rLJ">pic.twitter.com/PZVXqA0rLJ</a></p>&mdash; 神岡宇宙素粒子研究施設 (@Kamiokaobs_pr) <a href="https://twitter.com/Kamiokaobs_pr/status/1266215172403654659?ref_src=twsrc%5Etfw">May 29, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

## スーパーカミオカンデとは

[スーパーカミオカンデ](http://www-sk.icrr.u-tokyo.ac.jp/doc/sk/index.html)とは、岐阜県飛騨市の神岡鉱山内にある、ニュートリノや陽子崩壊の検出を目的とした装置、および実験の名前です。検出器は直径および高さが約 40 m の水タンクと、約10000本の光センサーからなっています。
詳しくは、[builderscon 2019](https://builderscon.io) において中の人である宇宙線研究所の早戸さんがエンジニア向けに話されている動画があるので、以下を参照してください。

- [スーパーカミオカンデの開発と運用 (Yoshinari Hayato) - builderscon tokyo 2019 - YouTube](https://www.youtube.com/watch?v=Fh8fIv8IE_0) 

## 公開されたデータについて

配布されている zip ファイルを展開すると以下のようになっており、CSV と PDF がそれぞれ8ファイル、Python スクリプトが1つ含まれています。CSV については[配布ページ](http://www-sk.icrr.u-tokyo.ac.jp/sk/detector/eventdisplay-data.html) に説明があります。 PDF については説明がないですが、おそらくそれぞれの CSV を図示したものと思われます。

```
1ring-e.0000.000002.csv
1ring-e.0001.000018.csv
1ring-mu.0000.000021.csv
1ring-mu.0001.000031.csv
2rings-pi0.0000.000065.csv
2rings-pi0.0001.000529.csv
multirings.0000.000014.csv
multirings.0001.000015.csv
run999999_001_000002.pdf
run999999_001_000014.pdf
run999999_001_000015.pdf
run999999_001_000018.pdf
run999999_001_000021.pdf
run999999_001_000031.pdf
run999999_001_000065.pdf
run999999_001_000529.pdf
skdata_sample.py
```

CSV はファイル名から4種類に分かれていることが分かりますが、名前から類推するに、`1ring-e` は電子ニュートリノの事象、`1ring-mu`はミューニュートリノの事象と思われます。`2rings-pi0` は、[こちらの粒子ID](https://particle.wiki/wiki/PDG_particle_numbering_scheme)などを参照して CSV を解読すると、ミューニュートリノが入ってきて中性のパイオンが生成され、パイオンが2つの光子に崩壊する事象と思われます。`multirings` については、ニュートリノが入ってきて複数の荷電粒子が生成された事象のようでした。[配布ページ](http://www-sk.icrr.u-tokyo.ac.jp/doc/sk/detector/eventdisplay-data.html)には「大気ニュートリノや陽子崩壊のシミュレーションデータを公開します」とあるのですが、陽子崩壊で3つのチェレンコフリングができる事象は含まれていないように見えます。今後に期待しましょう。（間違っていたらすみません）


## 図示してみる

### 先行例

2020/05/31 時点で以下の3つの先行例を観測しました。

付属している Python スクリプトを動かして図示する例は、すでに以下の記事で[ひだ宇宙科学館 カミオカラボ](https://www.city.hida.gifu.jp/site/kamiokalab/)の高知尾さんがやられているので、参照してください。

- [スーパーカミオカンデのシミュレーションデータで遊んでみた！｜ひだ宇宙科学館 カミオカラボ｜note](https://note.com/kamiokalab/n/n2c9ebfb8b0da)

また以下の記事では、名古屋大学の奥村さんが LEGO 風に描画されています。

- [Studio.app と LDraw を使って LEGO 風のイベントディスプレイをスーパーカミオカンデ用に作る - 宇宙線実験の覚え書き](https://oxon.hatenablog.com/entry/2020/05/31/122807)

KEK の高橋さんは2次元でプロットする方法を紹介されています。

- [Super Kamiokande のシミュレーションデータが公開されたので遊んでみた - Qiita](https://qiita.com/shotakaha/items/a86b98a39984b4e61b0a)

### アニメーションを作る

本記事では、 CSV に記録されている時間情報を使って、アニメーションを作ってみたいと思います。[Plotly の Python 版](https://plotly.com/python/)を使用して描いてみました。Plotly には [Plotly Express](https://plotly.com/python/plotly-express/) という簡単な API がありますが、そちらでは求めるアニメーションが作れなそうだったので、[plotly.graph_objects.Scatter3d](https://plotly.com/python/3d-scatter-plots/#basic-3d-scatter-plot) を使いました。以下、コードです。

```python
import pandas as pd
import sys
import plotly.graph_objects as go

print('input file :', sys.argv[1])
df = pd.read_csv(sys.argv[1], comment='#')
#df.columns = ['cable', 'Charge(p.e.)', 'time(ns)', 'x', 'y', 'z']
df.columns = ['cable', 'q', 't', 'x', 'y', 'z']

# 元のスクリプトと同じ規格化を行う
q_avg = df['q'].mean()
df['q_disp'] = [q_avg * 3 if q_val > q_avg * 3 else q_val for q_val in df['q']]
frame_interval = 50  # フレーム間を何 ns 空けるか。小さくすると重くなります。
frames = []
for i in range(int(min(df['t'])), int(max(df['t'])), frame_interval):
    # フレームの作成する。
    # i より前の時間のヒットは図に残しておく。
    df_slice = df[df["t"] < i]
    frames.append(
        go.Frame(data=go.Scatter3d(
            x=df_slice['x'], y=df_slice['y'], z=df_slice['z'])))

# アニメーションを作成する。
fig = go.Figure(
    data=[
        go.Scatter3d(
            x=df['x'],
            y=df['y'],
            z=df['z'],
            mode='markers',
            marker=dict(
                size=4,
                color=df['q_disp'],
                colorscale='Plotly3',  # カラースケールはお好みで
                #colorscale='Plasma',
                #colorscale='Bluered',
                opacity=1.0))
    ],
    layout=go.Layout(
        scene=dict(xaxis=dict(range=[min(df['x']), max(df['x'])],
                              autorange=False),
                   yaxis=dict(range=[min(df['y']), max(df['y'])],
                              autorange=False),
                   zaxis=dict(range=[min(df['z']), max(df['z'])],
                              autorange=False)),
        # autosize=True,
        title=sys.argv[1],
        updatemenus=[
            dict(type="buttons",
                 buttons=[dict(label="Play", method="animate", args=[None])])
        ]),
    frames=frames)

fig.write_html(sys.argv[1] + '.html', auto_open=False)  #アニメーションを HTML に書き込む。
# fig.show() # すぐに見たいときはコメントアウト
```

出来上がったアニメーションの例は以下です。

![mu21.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/63995/86bf685d-d984-0821-0d06-b2b6657b7a01.gif)

![multi14.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/63995/1fd1a797-d4c0-21e7-75f0-c9a225dbd5a5.gif)

![multi15.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/63995/805c6d71-5430-16da-c2ae-e55a53aa2a37.gif)


- ミューオンの事象 (1枚目) はリングがきれいに見えました。
- 軸のスケールがフレームごとに変わってしまうのは `go.Layout(scene=dict(xaxis= ...` の部分で指定しているはずなのですがうまくいきませんでした。

## Other References

- [Intro to Animations | Python | Plotly](https://plotly.com/python/animations/)
- [Built-in Continuous Color Scales | Python | Plotly](https://plotly.com/python/builtin-colorscales/)
- [Plotly.pyで3次元のムービーを作成 - Qiita](https://qiita.com/ufield/items/4c71809261e8addcf362)