Super Kamiokande のシミュレーションデータが公開されたので遊んでみた

# イベントディスプレイを作ってみよう！

- 5月29日（金）に、Super Kamiokande のFBページに下記の投稿がありました

> スーパーカミオカンデのイベントディスプレイ作成用データの公開を始めました。ニュートリノや陽子崩壊のシミュレーションデータを元に、あなた独自のイベントディスプレイを作ってみませんか？データには光電子増倍管の位置座標と光の量や時間の情報が含まれています。
> <del>http://www-sk.icrr.u-tokyo.ac.jp/sk/detector/eventdisplay-data.html</del>

（2024/02/22追記）
ひさびさに確認したらURLが変更されていました。
現在は以下のURLページの下の方に説明があります。

https://www-sk.icrr.u-tokyo.ac.jp/sk/about/detector/

# サンプルデータをダウンロードする

- データは下記の 東京大学宇宙線研究所 神岡宇宙素粒子研究施設 のウェブページにある
  - http://www-sk.icrr.u-tokyo.ac.jp/sk/detector/eventdisplay-data.html
- 「サンプルデータのダウンロードはこちらから」をクリックしてダウンロード
- ``zip``ファイルになっているので、展開して適当な場所に移動させる
  - 以下では、``~/Downloads/sample/`` にあるつもりで進める

## ダウンロードしたファイルの概要

- ``CSV``ファイル : お絵かきに使うためのファイル一式
  - CSVファイルのセルの内容についてはSKのウェブサイトを参照
  - ファイル名は、どのようなシミュレーションをしたのか、少しわかるようになっている
    - ``1ring-e.*.csv`` : **電子**が作ったチェレンコフリングが**1つだけ**
    - ``1ring-mu.*.csv`` : **ミューオン**が作ったチェレンコフリングが**1つだけ**
    - ``2rings-pi0.*.csv`` : $\pi^0$ 中間子が作ったチェレンコフリングが**2つ**
    - ``multirings*.csv`` : チェレンコフリングが**たくさん**
- ``PDF``ファイル
  - ``CSV``ファイルに対応したイベントディスプレイ
  - 自分でイベントディスプレイを作るときの参考にしたり、答え合わせに使ったらいいと思う
- ``Python``ファイル
  - ``skdata_samply.py`` : SKの研究者が作ったサンプルプログラム
  - 特に説明はないけれど、中身を読めば何をしているか、だいたい分かるように書かれている

# サンプルプログラムを実行

- とりあえず動かそう！
- 以下のように ``CSVファイル名`` を引数に指定して、サンプルプログラムを実行する

```bash
$ python skdata_sample.py CSVファイル名
```

## 実際に入力してみたコマンド

- マウスでぐりぐりできる3Dプロットがでてくる

```bash
$ cd ~/Downloads/sample/
$ python3 skdata_sample.py 1ring-e.0000.000002.csv
```

---

ここまでで、サンプルプログラムの使い方を簡単に確認できたので、遊んでみた

# 遊んでみた

- ``pandas`` + ``altair`` で作ってみた
- ``JupyterLab``に以下のコマンドを入力していけばOK

## 実行した環境

- ``Python`` : ``3.7.3``
- ``JupyterLab`` : ``2.1.3``
- ``Pandas`` : ``1.0.4``
- ``Altair`` : ``4.1.0``

## 編集＆追記（2020/06/01）

- 大気ニュートリノのデータがたくさん（合計300イベント）追加された
- 雑に書いていたので、この記事をいろいろ編集
  - コードを書き直したり
  - プロットを新規に追加したり
  - 物理の説明を追記したり
- GitLab にリポジトリを作成
  - https://gitlab.com/shotakaha/skdisplay
  - ``skdisplay.ipynb`` を参照
  - ``README.md`` では物理も説明している


# ファイル名を取得する

- ダウンロードしたデータは ``skdata`` 以下の保存することにした
- ファイル名をいちいち確認して指定するのはめんどくさいので、``pathlib.Path`` モジュールを使って ``*.csv`` ファイル名の一覧を ``リスト型`` で作ることにした

## ``pathlib``モジュール


- ``pathlib.Path`` : https://docs.python.org/3/library/pathlib.html
  - ``pathlib.Path`` モジュールを使いたいので ``from pathlib import Path`` する
- ``pathlib.Path.glob`` : https://docs.python.org/3/library/pathlib.html#pathlib.Path.glob
  - サブディレクトリの中にあるファイルを取得するには ``.glob('*/*.csv')`` とする必要がある
- ``pathlib.Path.rglob`` : https://docs.python.org/3/library/pathlib.html#pathlib.Path.rglob
  - サブディレクトリの中にあるファイルを取得するのも ``.rglob('*.csv')`` だけでOK


```python
from pathlib import Path
datad = Path('skdata')
fnames = sorted(datad.glob('*/*.csv'))
#fnames = sorted(datad.rglob('*.csv'))
```

- 読み込んだファイルの数を確認

```python
len(fnames)
```

# ``Pandas`` を使ってデータフレームとして読み込む

- まずCSVファイルの中身を確認する
  - 最初の数行はヘッダ情報で # から始まっている（# は C++だとコメントに使う文字）
  - 列のラベルは ``cable``（ケーブル番号）、``charge``（電荷[p.e.]）、``time``（時間[ns]）、``x``（[mm]）、``y``（[mm）、``z``（[mm]）
- CSVの読み込みには pandas.read_csvを使う
  - https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.read_csv.html
  - ヘッダ情報はファイルによって行の数が異なるので ``skiprows`` オプションは使えない
  - 他にオプションを探してみると ``comment`` オプションが使える（今回初めて知った）
  - 列のラベルも自分で指定する必要があるので names オプションも使う
- とりあえずこれくらいでよさそう
- ``ファイル名`` を引数にして、データフレームを作成する関数を作った

```python
import pandas as pd

def load(fname):
    names = ['cable', 'charge', 'time', 'x', 'y', 'z']
    data = pd.read_csv(fname, comment='#', names=names)
    print(f'Loaded! {fname}')
    return data
```

```python
fname = files[0]  ## 読み込むファイルを適当に指定
data = load(fname)
# ==> Loaded! skdata/1ring-e.100_events/1ring-e.0000.000002.csv
```

# ``Altair`` を使って 2Dの断面図 を作成する

- ``Altair`` : https://altair-viz.github.io/
- ``Example Gallery`` から今回の用途に合いそうなサンプルを探して参考にする
  - https://altair-viz.github.io/gallery/index.html
- 断面図を作成するには、散布図を使えばよさそう
  - https://altair-viz.github.io/gallery/binned_scatterplot.html
- プロットにマウスオーバーすると、そのマーカーの情報を表示することができる
  - https://altair-viz.github.io/gallery/scatter_tooltips.html

```python
import altair as alt

def crosssection(data, x, y, color, title, w=300, h=300):
    tips = list(data.columns)
    chart = alt.Chart(data).mark_circle().encode(
        x=x,
        y=y,
        color=color,
        tooltip=tips,
    ).properties(
        title=title,
        width=w,
        height=h,
    )
    return chart
```

- ``x-y`` の断面図
- ``y-z`` の断面図
- ``x-z`` の断面図
- ``color`` は ``電荷``

```python
color = 'charge'
xy = crosssection(data, x='x', y='y', color=color, title='X-Y 断面')
yz = crosssection(data, x='y', y='z', color=color, title='Y-Z 断面')
xz = crosssection(data, x='x', y='z', color=color, title='X-Z 断面')
```

- 横に並べてプロットを表示

```python
xy | yz | xz
```

![skdisplay01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/273122/11eff733-722e-e595-b664-b8316fe77f63.png)

- （真ん中の図で）検出器の右下のほうが光ってるのが分かる
- ``VSCode``／``Jupyter Lab`` などでプロットを作っていると、マウスをホバーさせたときに、そのポイントの情報が表示される（ここに貼ってあるのはただのPNGなので何もおきない）

## 見やすくするために改善した（2020/06/07追記）

- 検出器の形状を考えると、この表示方法では、**手前のヒット**と**奥のヒット** が重なってしまっている
- 見やすくするために、``Z軸`` に使っているラベルの値で手前／奥を選択して、ディスプレイを描きわけることにする
- 描きわける方法は、大きく分けて以下のタイミングでできる
  - ``pandas`` のデータフレームで選別する
  - ``altair`` でプロットを表示するときに選別する
- 個人的には ``pandas`` でやるほうがよい（と思っている）
  - ``pandas`` を使ったやり方も何種類か思うが、今回は ``pandas.DataFrame.query`` を使う
  - https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.query.html
- 全部で6種類のデータを用意する
  - queried data という意味で、変数名は ``qdata`` とした
  - プロットを表示する関数（``crosssection``）はそのまま使える

```python
color = 'charge'

# 1.
qdata = data.query('z >= 0')    
xy1 = crosssection(qdata, x='x', y='y', color=color, title='X-Y 断面（Z>0）')
# 2.
qdata = data.query('z < 0')    
xy2 = crosssection(qdata, x='x', y='y', color=color, title='X-Y 断面（Z<0）')
# 3.
qdata = data.query('x >= 0')    
yz1 = crosssection(qdata, x='y', y='z', color=color, title='Y-Z 断面（X>0）')
# 4.
qdata = data.query('x < 0')    
yz2 = crosssection(qdata, x='y', y='z', color=color, title='Y-Z 断面（X<0）')
# 5.
qdata = data.query('y >= 0')        
xz1 = crosssection(qdata, x='x', y='z', color=color, title='X-Z 断面（Y>0）')
# 6.
qdata = data.query('y < 0')
xz2 = crosssection(qdata, x='x', y='z', color=color, title='X-Z 断面（Y<0）')
```

- 6つのプロットを 3列×2行にならべた
  - ``|`` で列に並べることができる
  - ``&`` で行に並べることができる

```python
# ``Z軸の値 > 0`` のプロットを上段に並べる
# ``Z軸の値 < 0`` のプロットを下段に並べる
(xy1 | yz1 | xz1) & (xy2 | yz2 | xz2)
```

![skdisplay02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/273122/aa9b6499-a9bd-6ba9-aac3-aa5b76c48ba6.png)

- 直前のプロットと比べると、少し見やすくなった


# ここまでをまとめる

- ファイルの読み込み : ``load(fname)``
- 断面図の作成 : ``crosssection(data, x, y, color, title)``

```python
def display(fname, color='charge'):
    data = load(fname)
    xy = crosssection(data, x='x', y='y', color=color, title='X-Y 断面')
    yz = crosssection(data, x='y', y='z', color=color, title='Y-Z 断面')
    xz = crosssection(data, x='x', y='z', color=color, title='X-Z 断面')
    chart = xy | yz | xz
    return chart
```

- ファイル名を適当に与えれば、プロットが表示される（はず）

---

- **電子が作った**チェレンコフリングが **1つだけ** のシミュレーションデータのプロット例

```python
display(fnames[20])
# ==> Loaded! skdata/1ring-e.100_events/1ring-e.0020.000407.csv
```

![skdisplay020.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/273122/5cb05206-4509-9014-247f-0afa9baa1f05.png)

---

- **ミューオンが作った**チェレンコフリングが **1つだけ** のシミュレーションデータのプロット例

```python
display(fnames[146])
# ==> Loaded! skdata/1ring-mu.100_events/1ring-mu.0045.000496.csv
```

![skdisplay146.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/273122/be8ea154-7ec3-e075-a891-e42a36ff248e.png)

---

- チェレンコフリングが **複数ある** シミュレーションデータのプロット例

```python
display(fnames[233])
# ==> Loaded! skdata/multirings.100_events/multirings.0031.000910.csv
```

![skdisplay233.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/273122/db91dbe4-2a4d-d581-8061-9957c302d697.png)


# まとめ

たのしい！😊

---

> 物理のことも少しずつ書き足していきます

# スーパーカミオカンデについて

## スーパーカミオカンデの形状

- 直径：約40m、高さ：約40m の円柱型の検出器
- 検出器の内部は全部『水』で満たされている（5万トンの超純水を使用）
- 検出器の内側には超高感度の光検出器が約13000本取り付けてある
- 今回配布されているデータの座標 (x, y, z) は、この光検出器が取り付けている場所を示している
- 詳細は [検出器の構造 - SK](http://www-sk.icrr.u-tokyo.ac.jp/sk/detector/) を参照

## スーパーカミオカンデで検出しているもの

- 検出器の内側に取り付けた光検出器で、素粒子反応から生じる **チェレンコフ光** を検出
- 光検出器は **光** を **光電子** に変換＆増幅し、**電気信号** に変換する装置
- 今回配布されているデータの
  - 電荷（``charge (p.e.)``）は、光電子の数（単位 ``p.e.`` は photo-electron（光電子）の略）
  - 時間（``time (ns)``） はイベントが発生してから、光検出器にチェレンコフ光が到達するのにかかった時間（ナノ秒）
- 詳細は [検出方法 - SK](http://www-sk.icrr.u-tokyo.ac.jp/sk/detector/cherenkov.html) を参照


## チェレンコフ光

- 真空中の **光速c** は 約30万 km/s
- **物質中の光速** は **屈折率に比例** する（ **1／屈折率** になる＝遅くなる）
- 真空中の光速を超えることはできない（＝光速度不変の原理）が、**物質中の光速** は超えることができる
- 粒子が**物質中の光速**を超えるとき、進行方向に**円錐状に電磁波**を発する
- これを**チェレンコフ光**と呼ぶ
  - チェレンコフさん（7月28日生まれ）が発見
- 水の場合、屈折率が **約1.3** なので、水中の光速は （30万/1.3）= **約23万 km/s** と計算できる
- アクリルの場合、屈折率が **約1.5** なので、アクリル中の光速は （30万/1.5）= **約20万 km/s** と計算できる


## ニュートリノの種類と反応

- スーパーカミオカンデは素粒子『ニュートリノ』を研究する大型の実験装置
- ニュートリノは3種類あることが分かっている
  - **電子型ニュートリノ** : 反応後に電子が飛び出てくる
  - **ミュー型ニュートリノ** : 反応後にミューオンが飛び出てくる
  - **タウ型ニュートリノ** : 反応後にタウが飛び出てくる
- また、作られる方法によっても呼び方が分けられている
  - **太陽ニュートリノ** : 太陽の核融合反応で作られて、地球まで飛んでくるニュートリノ
  - **大気ニュートリノ** : 1次宇線が地球の大気がぶつかって生じるニュートリノ
  - **加速器ニュートリノ** : 大型の加速器を使って人工的に作るニュートリノ
  - （他にもある）