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

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

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

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

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

データは下記の東京大学宇宙線研究所神岡宇宙素粒子研究施設のウェブページにある
- 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ファイル名 を引数に指定して、サンプルプログラムを実行する

$ python skdata_sample.py CSVファイル名

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

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

$ 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

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

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

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 オプションも使う
とりあえずこれくらいでよさそう
ファイル名 を引数にして、データフレームを作成する関数を作った

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

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

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 は 電荷

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 断面')

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

xy | yz | xz

（真ん中の図で）検出器の右下のほうが光ってるのが分かる
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）はそのまま使える

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行にならべた
- | で列に並べることができる
- & で行に並べることができる

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

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

ここまでをまとめる

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

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つだけ のシミュレーションデータのプロット例

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

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

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

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

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

まとめ

たのしい！😊

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

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

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

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

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

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

チェレンコフ光

真空中の 光速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次宇線が地球の大気がぶつかって生じるニュートリノ
- 加速器ニュートリノ : 大型の加速器を使って人工的に作るニュートリノ
- （他にもある）

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

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

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

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

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

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

遊んでみた

実行した環境

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

ファイル名を取得する

pathlibモジュール

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

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

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

ここまでをまとめる

まとめ

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

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

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

チェレンコフ光

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

`pathlib`モジュール

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

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