Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
0
Help us understand the problem. What is going on with this article?
@shotakaha

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

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

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

スーパーカミオカンデのイベントディスプレイ作成用データの公開を始めました。ニュートリノや陽子崩壊のシミュレーションデータを元に、あなた独自のイベントディスプレイを作ってみませんか?データには光電子増倍管の位置座標と光の量や時間の情報が含まれています。
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 にリポジトリを作成

ファイル名を取得する

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

pathlibモジュール

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

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

skdisplay01.png

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

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

  • 検出器の形状を考えると、この表示方法では、手前のヒット奥のヒット が重なってしまっている
  • 見やすくするために、Z軸 に使っているラベルの値で手前/奥を選択して、ディスプレイを描きわけることにする
  • 描きわける方法は、大きく分けて以下のタイミングでできる
    • pandas のデータフレームで選別する
    • altair でプロットを表示するときに選別する
  • 個人的には pandas でやるほうがよい(と思っている)
  • 全部で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)

skdisplay02.png

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

ここまでをまとめる

  • ファイルの読み込み : 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

skdisplay020.png


  • ミューオンが作ったチェレンコフリングが 1つだけ のシミュレーションデータのプロット例
display(fnames[146])
# ==> Loaded! skdata/1ring-mu.100_events/1ring-mu.0045.000496.csv

skdisplay146.png


  • チェレンコフリングが 複数ある シミュレーションデータのプロット例
display(fnames[233])
# ==> Loaded! skdata/multirings.100_events/multirings.0031.000910.csv

skdisplay233.png

まとめ

たのしい!😊


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

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

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

  • 直径:約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次宇線が地球の大気がぶつかって生じるニュートリノ
    • 加速器ニュートリノ : 大型の加速器を使って人工的に作るニュートリノ
    • (他にもある)
0
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
shotakaha
素粒子物理学実験してました。とある研究機構の科学広報(仮)です。合気道、自転車、カメラ、登山などします。広く浅く。みでしみん。つくば科学教育マイスター。素粒子のお兄さん。あっぷちゃん。野望:サイエンスゼ〇に出ること。

Comments

No comments
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account Login
0
Help us understand the problem. What is going on with this article?