物理の駅 by onsanai

Physics Station → PhSt 質問・疑問・間違いの指摘は、コメントに書くか、直接伝えるときっと良いことがあります。主にWindows or Ubuntu用の記事です

単スリットのフラウンホーファー回折のシミュレーション

http://www.sci.keio.ac.jp/gp/87B7D75A/A6070F75/6E345155.pdf

単スリットについて考えてみた。 スリット幅が波長に比べて十分広いか、スクリーンの位置が十分遠い場合には、回折による干渉模様を観察できる。

数式をPDFファイルからパクって書いておく。

\displaystyle
\begin{eqnarray}
\frac{xd'}{l\lambda}
 =
  \begin{cases}
    0, \pm\frac{3}{2}, \pm\frac{5}{2}, \pm\frac{7}{2},\cdot\cdot\cdot のとき明線 \\
\\
    \pm 1, \pm 2, \pm 3 \cdot\cdot\cdot  のとき暗線
  \end{cases}
\end{eqnarray}


  • x: スクリーン中央を0としたときの各干渉模様の位置
  • d: スリット間隔 今回は0.01 mm、0.1 mm、1 mmを試してみる
  • l: スリットからスクリーンまでの距離 干渉模様の位置を同じにするため、10 mm、100 mm、1000 mmを試してみる
  • λ: レーザー光の波長 そんなのないけど500 nmの波長とする

シミュレーションではスクリーンの空間で2ミクロンごと(画像の1ピクセルが2ミクロン相当、縦横2048ピクセル)に計算させた。最小のスリット間隔 10ミクロンの1/5なので十分だと思う(たぶん)。

f:id:onsanai:20200125031334j:plainf:id:onsanai:20200125031337j:plainf:id:onsanai:20200125031351j:plain
左から順にスリット間隔0.01mm、0.1mm、1mm

フラウンホーファー回折が成り立つ条件は、ざっくり

\displaystyle
\frac{d^2}{l\lambda}\ll1

らしい。各条件での左辺の値は0.02、0.2、2、であり、確かに一番右の画像ではきれいな回折像が確認できなくなっている。同じスリット幅でも、スリットとスクリーンまでの距離を10倍に離すと、以下の像が得られる。

f:id:onsanai:20200125033639j:plain:w220

ただし、横方向の干渉模様の位置を同じにするため、画像の1ピクセルを20ミクロン相当にした。そのため、縦方向にも縮んでしまっているが、横方向の回折像はくっきり確認できる。

学生による実験レポートのコピペについて

実験レポートの(やってはいけない)コピペは、起きてほしくはないが、よく起こる。

レポートのクオリティを上げたいとき、比較的よく書けている他人の実験レポートの一部や全部をコピペし、間違いを修正し、加筆するのは、クオリティを上げるという観点だけで言えば合理的な手段である。また、全く同じ実験データを用いて結果を示す時、表やグラフを再利用することは、やはり合理的な手段である。

自分の研究を振り返ってみると、複数の研究者による共同実験で、同じ時期に異なる人が外部に対して発表する時に、全ての発表資料を一から作る必要があるかというとそんなことはない。比較的よく作られた発表資料をコピペし、自分の発表用に加筆修正することは、合理的で妥当な手段である。全く同じ実験データを用いて異なる論文を書くとき、論文ごとに表やグラフを再利用することは、著作権の問題さえクリアすればやはり合理的で妥当な手段である。

学生実験は、2~3人で実験を行うことが多く、その結果は同じになることが多い。 よって、複数の研究者による共同実験でありがちな、表やグラフのコピペは、節度とルールを守れば、許容されると思う。

ここで大事なのは、レポートの地の文(表やグラフ以外の本分)である。地の文は、独立に書けば同じになることはまずない。導入、実験手法、結果、考察に至るまで、同じ結果を用いたとしても、同じ文章になることは絶対にない。よって、地の文で他人の文章を真似ることに対しては、厳しく処罰するべきだろう。以降、この地の文で他人の文章を真似する行為をやってはいけないコピペ(以下、コピペ)と定義することにする。

実験レポートの目的に立ち戻れば、学生にとってレポートの質を上げることは、レポート点(?)を上げることに繋がり、学生の直接的な利益になる。すぐに効果が実感できないが、レポートを書く技術の向上は、現在と近い将来の学業・仕事においても利益となる。教える側にとって、学生がコピペすることなくレポートの質が上がることは、教える側の本分である。

それらを踏まえても、やはりコピペは発生するものである。なので、色んなコピペ対策を検討してみた。私の感想なので、他の見解がある場合は一緒に議論させて欲しい。

コピペをさせないために手書きでレポートを提出させる

一見効果がありそうだが、実はあまり意味がない。なぜなら手書きでもコピペは出来るからだ。学生側に実験レポートの蓄積があり、そのいずれかからコピペされた場合は、読み手にコピペを確認する術はない。コピペはどうにもならないので、手書きでもいいから一度は書かせたほうがまし、みたいな大学があるのは知っているが、それはどうなんだろう。あと、読む側にとって読み辛い手書きも多々ある。

コピペチェッカーを導入する

高価なので却下、と言いたいところだが、今どきは頑張れば作れちゃう気もする。だが、そこで頑張らなくても良いし、頑張りたくない。

コピペはやるなと言う

やらないことは当然のことなので、言っても意味がない気がする。規範意識とはその人の持つ素質なので、言ってもやる奴には何を言っても無駄な気がする。

コピペはバレることを教える

バレるからコピペをするしないの話ではないので、目的と手段が一致していないので良くない。また、バレないように過去のレポートからコピペしだすととても面倒である。

連帯責任

コピペ元や、学年全体に罰を与えることは、相互信頼を醸成すべき学生生活において、意味がないばかりか害である。そもそも、コピペは本人の問題だし、学生のコピペを許してしまうことは担当の教員の責任なので、他の学生にその責任を負わせることは完全に筋違いである。

コピペをすると落第させることを通知する

実験の単位を認めないことはほぼ確実に留年(=落第)を意味するため、かなり強い罰になり、効果は期待できる。ただ、落第させると通知しても、バレないようにコピペをしようとする学生はいるかもしれない。

何がコピペであるかを教える

これはとても重要だ。新入の学生にとって、コピペの定義は教わらないと分からない。

レポートを書く上で楽をする正しいやり方を教える

これも有効かもしれない。共同実験者の表やグラフを、適切に再利用することが許されることは教えたほうが良いだろう。

なぜレポートを書くのかについて教える、学生間で議論する

自分の行為の目的を教えることや考えることは重要だ。ただし、目的を理解することが行為に結びつくかどうかは本人次第でもある。

レポート内容、特に考察について口頭試問を行う

口頭試問を行うと、レポートがコピペであるかどうかを簡単に確認できる。ただし、口頭試問を行う側の負担は大きく、能力を必要とし、そしてめちゃくちゃ面倒くさい。

ダブルΛハイパー核 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

f:id:onsanai:20191021023355p:plain
Weightあり

f:id:onsanai:20191021023338p:plain
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

直交座標系のずれ量を飛跡の角度に相当するRadial Lateral空間に変換するコード

放射線の飛跡のように、座標(cx, cy) 角度 (ax, ay)で表すことの出来る量の比較をする時、しばしば飛跡の進行方向(Radial方向)と垂直方向(Lateral方向)に分離して考えることがある。

飛跡1 (cx1, cy1, ax1, ay1)、飛跡2 (cx2, cy2, ax2, ay2) があるとき、角度差は次のように求めることができる。

import math

# 方法1
def xy2radial(cx1,cy1,ax1,ay1,cx2,cy2,ax2,ay2):
    ax = (ax1 + ax2) * 0.5
    ay = (ay1 + ay2) * 0.5
    r = (ax ** 2 + ay ** 2) ** 0.5
    ax /= r
    ay /= r
    dax = ax2 - ax1
    day = ay2 - ay1
    dar = ax * dax + ay * day
    dal = -ay * dax + ax * day
    dcx = cx2 - cx1
    dcy = cy2 - cy1
    dcr = ax * dcx + ay * dcy
    dcl = -ay * dcx + ax * dcy

    return dcr,dcl,dar,dal

# 方法2
def xy2radial2(cx1,cy1,ax1,ay1,cx2,cy2,ax2,ay2):
    ax = (ax1 + ax2) * 0.5
    ay = (ay1 + ay2) * 0.5
    rotation_cos = math.cos(-math.atan2(ay,ax))
    rotation_sin = math.sin(-math.atan2(ay,ax))
    dax = ax2 - ax1
    day = ay2 - ay1
    dar = dax * rotation_cos - day * rotation_sin
    dal = dax * rotation_sin + day * rotation_cos
    dx = cx2 - cx1
    dy = cy2 - cy1
    dcr = dx * rotation_cos - dy * rotation_sin
    dcl = dx * rotation_sin + dy * rotation_cos

    return dcr,dcl,dar,dal

処理時間は若干方法1の方が早かった。三角関数の処理は、CPUで実装されている? らしいのでそんなに遅くもない。

f:id:onsanai:20190808123417p:plain

参考文献 https://dx.doi.org/10.1093/ptep/ptx131