物理の駅 Physics station by 現役研究者

テクノロジーは共有されてこそ栄える

Pytyon: Jupyter lab 起動時に cannot import name 'soft_unicode' from 'markupsafe'

Python Jupyter 
C:\Users\Masahiro\Downloads\slack>jupyter lab
Traceback (most recent call last):
  File "c:\users\masahiro\appdata\local\programs\python\python37\lib\runpy.py", line 193, in _run_module_as_main
    "__main__", mod_spec)
  File "c:\users\masahiro\appdata\local\programs\python\python37\lib\runpy.py", line 85, in _run_code
    exec(code, run_globals)
  File "C:\Users\Masahiro\AppData\Local\Programs\Python\Python37\Scripts\jupyter-lab.EXE\__main__.py", line 4, in <module>
  File "c:\users\masahiro\appdata\local\programs\python\python37\lib\site-packages\jupyterlab\__init__.py", line 7, in <module>
    from .labapp import LabApp
  File "c:\users\masahiro\appdata\local\programs\python\python37\lib\site-packages\jupyterlab\labapp.py", line 12, in <module>
    from jinja2 import Environment, FileSystemLoader
  File "c:\users\masahiro\appdata\local\programs\python\python37\lib\site-packages\jinja2\__init__.py", line 12, in <module>
    from .environment import Environment
  File "c:\users\masahiro\appdata\local\programs\python\python37\lib\site-packages\jinja2\environment.py", line 25, in <module>
    from .defaults import BLOCK_END_STRING
  File "c:\users\masahiro\appdata\local\programs\python\python37\lib\site-packages\jinja2\defaults.py", line 3, in <module>
    from .filters import FILTERS as DEFAULT_FILTERS  # noqa: F401
  File "c:\users\masahiro\appdata\local\programs\python\python37\lib\site-packages\jinja2\filters.py", line 13, in <module>
    from markupsafe import soft_unicode
ImportError: cannot import name 'soft_unicode' from 'markupsafe' (c:\users\masahiro\appdata\local\programs\python\python37\lib\site-packages\markupsafe\__init__.py)

qiita.com

を参考に、pip install Flask で解決した。

Python: 2次元図で凸多角形の内部のイベントを得る

Python

この記事の Awkward のバージョンは1.X.Xです

例外処理はちゃんとしてないので注意せよ。

import numpy as np
import awkward as ak
import matplotlib.pyplot as plt

# 配列作成
obj={}
rs = np.random.RandomState(0)
obj["memx"]=rs.rand(1000)
obj["memy"]=rs.rand(1000)
memx = "memx"
memy = "memy"
tree = ak.Array(obj)

# グラフ化
plt.scatter(tree[memx],tree[memy])
plt.plot(pointx,pointy,color="tab:red")
plt.xlim(-0.1,1.1)
plt.ylim(-0.1,1.1)
plt.show()

def get_in_convex_polygon(tree, memx, memy, pointx, pointy):
    assert(len(pointx)==len(pointy))
    assert(pointx[0]==pointx[-1])
    assert(pointy[0]==pointy[-1])
    assert(len(pointx)>=4)
    
    dx = pointx[1]-pointx[0]
    dy = pointy[1]-pointy[0]
    if ((pointx[2]-pointx[0])*dy-(pointy[2]-pointy[0])*dx)>0:
        # 回転方向を合わせる
        pointx.reverse()
        pointy.reverse()
    
    for i in range(len(pointx)-1):
        dx = pointx[i+1]-pointx[i+0]
        dy = pointy[i+1]-pointy[i+0]
        tree = tree[((tree[memx]-pointx[i])*dy-(tree[memy]-pointy[i])*dx)<0]
    return tree

# 三角形の点列 (終点は始点と同じにせよ)
pointx = [0.0,0.6,0.3,0.0]
pointy = [0.0,0.0,0.5,0.0]

tree = get_in_convex_polygon(tree,memx,memy,pointx,pointy)

# グラフ化
plt.scatter(tree[memx],tree[memy])
plt.plot(pointx,pointy,color="tab:red")
plt.xlim(-0.1,1.1)
plt.ylim(-0.1,1.1)
plt.show()

凸多角形って言うけど、凸って字は凸多角形ではないんだよなあ。

Ubuntu 18.04 LTS 起動しない問題の解決 (ディスク関連)

Linux/Ubuntu

「Ubuntuが起動しないので何とかしてほしい」という依頼を受けて、OSが入ったSystem用のSSDだけもらって原因を究明した話。

手元のデスクトップパソコンに適当に挿して起動すると

/dev/sda2: clean, **/** files, **/** blocks
You are in emergency mode. After logging in, type "journalctl -xb" to view
system logs, "systemctl reboot" to reboot, "systemctl default" or "exit" 
to boot into default mode.
Press Enter for maintenance
(or press Control-D to continue):

と出た。Enterでメンテナンスモードに入る。

journalctl -xb で、起動時のシステムログが見れるので、とりあえず見ていく。赤字のところがエラーっぽいところ。

Disk1がマウントできない? そういえば、マウントしてたなあと思い出して、起動時にマウントするファイルシステムの情報を見る。コマンドは sudo pico /etc/fstab。

Disk1がついてた。今は必要ないので、この行の先頭に#をつけてコメントアウトする。

念のため、ディスクの健全性も確認する。8GB以上の容量のデータを消しても良いUSBフラッシュメモリを用意する。

Ubuntu 18.04 LTSのisoをダウンロードして、

www.ubuntulinux.jp

RufusでUSBフラッシュメモリに書き込む。

rufus.ie

これを刺して、BIOSのBootでUSBフラッシュメモリを選び起動すると Try Ubuntu without installing が選べるのでそれを選ぶ。

トライアル版のUbuntuが立ち上がるので、Disk(ディスク)というアプリを立ち上げ、Ubuntuがインストールされているディスクを選び、ボリュームを選び、右クリックで「Check Filesystem」とする。

EFIシステムの方で「ファイルシステムは損傷しています 修復が必要です」というエラーが出たので、「Repair Filesystem」とする。

「ファイルシステム修復時にエラーが発生しました」という怖いメッセージが出たが、どうも修復はできた模様。(なぜならば、もう一度起動してCheck Filesystemするとエラーが出なかったので)

このまま、トライアル版のUbuntuはシャットダウンして、USBフラッシュメモリを抜いて起動すると、当方の環境だと起動できました。

なお、実際の作業では、Windows PCにSystem用SSDを挿して、CrystalDiskInfoでSMART情報を確認してSSDが壊れていないことを確認→トライアル版のUbuntuのDiskアプリでファイルシステムを修復→マウント情報の整理→起動 という流れでした。

Python: 重イオンのEnergy loss と Energy loss straggling

Python 数学/物理

まず、各定数の確認する。

import numpy as np

elementary_charge = 1.602176634*10**-19 # 電気素量 C
c_light_m = 299792458 # 光速 m/s
electron_mass_kg = 9.1093837015*10**-31 # 電子質量 kg
permittivity_vacuum = 8.8541878128*10**-12 # 真空の誘電率 F/m = C/Vm
r_e = 2.8179403262*10**-13 #古典電子半径 cm
print("r_e =", elementary_charge**2 / (4 * np.pi * permittivity_vacuum * electron_mass_kg * c_light_m**2)*100,"cm")

electron_mass = 0.510998928 #   // 電子質量 MeV/c^2
Na = 6.02214179*10**23 # アボガドロ数 /mol
dedx_constant = 0.3070749187 #4*pi*Na*me*c^2*r_e^2  //MeV cm^2
print("dedx_constant =", 4 * np.pi * Na * electron_mass * r_e**2,"MeV/cm2")

色んな項を無視したエネルギー損失と、エネルギーストラグリングを計算する。

def energy_loss_straggling(p_A, p_Z, p_Q, p_T, t_Z, shell_correction = False, LSX=1):
    '''
    p_A=238 # projectileのA
    p_Z=92 # projectileのZ
    p_Q=92 # projectileのQ
    p_T = 345 #projectileのEnergy MeV/u
    t_Z = 4       # tagetのZ
    '''
    import numpy as np

    electron_mass = 0.510998928 #   // 電子質量 MeV/c^2
    Na = 6.02214179*10**23 # アボガドロ数 /mol
    dedx_constant = 0.3070749187 #4*pi*Na*me*c^2*r_e^2  //MeV cm^2
    atomic_mass_unit = 931.4940954 # // 統一原子質量単位 MeV/c^2

    zp_eff = p_Q #projectileのQ
    gamma=1.0 + p_T/atomic_mass_unit
    beta2=1.0-1.0/(gamma*gamma)
    beta = np.sqrt(beta2)

    def get_a(t_Z):
        standard_weights = [1.008, 4.003,\
            6.941, 9.012, 10.81, 12.01, 14.01, 16.00, 19.00, 20.18,\
            22.99, 24.31, 26.98, 28.09, 30.97, 32.07, 35.45, 39.95,\
            39.10, 40.08, 44.96, 47.87, 50.94, 52.00, 54.94, 55.85, 58.93, 58.69, 63.55, 65.38, 69.72, 72.63, 74.92, 78.97, 79.90, 83.80,\
            85.47, 87.62, 88.91, 91.22, 92.91, 95.95, 99, 101.1, 102.9, 106.4, 107.9, 112.4, 114.8, 118.7, 121.8, 127.6, 126.9, 131.3,\
            132.9, 137.3,\
            138.9, 140.1, 140.9, 144.2, 145, 150.4, 152.0, 157.3, 158.9, 162.5, 164.9, 167.3, 168.9, 173.0, 175.0,\
            178.5, 180.9, 183.8, 186.2, 190.2, 192.2, 195.1, 197.0, 200.6, 204.4, 207.2, 209.0, 210, 210, 222,\
            223, 226,\
            227, 232.0, 231.0, 238.0, 237, 239, 243, 247, 247, 252, 252, 257, 258, 259, 262,\
            267, 268, 271, 272, 277, 276, 281, 280, 285, 278, 289, 289, 293, 293, 294]
        return standard_weights[t_Z-1]

    t_A = get_a(t_Z)

    def get_ipot(t_Z):
        arr_ipot = [
            19.2, 41.8, \
            40.0, 63.7, 76.0, 81.0, 82.0, 95.0, 115.0,137.0,\
            149.0,156.0,166.0,173.0,173.0,180.0,174.0,188.0,\
            190.0,191.0,216.0,233.0,245.0,257.0,272.0,286.0,297.0,311.0,322.0,330.0,334.0,350.0,347.0,348.0,343.0,352.0,\
            363.0,366.0,379.0,393.0,417.0,424.0,428.0,441.0,449.0,470.0,470.0,469.0,488.0,488.0,487.0,485.0,491.0,482.0,\
            488.0,491.0,501.0,523.0,535.0,546.0,560.0,574.0,580.0,591.0,614.0,628.0,650.0,658.0,674.0,684.0,694.0,705.0,718.0,727.0,736.0,746.0,757.0,790.0,790.0,800.0,810.0,823.0,823.0,830.0,825.0,794.0,\
            827.0,826.0,841.0,847.0,878.0,890.0,900.0,910.0,920.0,930.0,940.0,950.0,960.0,970.0,980.0,990.0,1000.,1010.,1020.,1030.,1040.,1050.,1060.,1070.,1080.,1090.,1100.,1110.,1120.,1130.,1140.,1150.
            ]
        return arr_ipot[t_Z-1]
    Ipot = get_ipot(t_Z) #平均イオン化ポテンシャル [eV]
    
    f1 = dedx_constant*pow(zp_eff,2.0)*t_Z/(beta2*t_A)
    f2 = np.log(2.0*electron_mass*beta2*gamma*gamma/(Ipot/1000000)) - beta2
    
    cor = 0 #shell補正
    f2 = f2 - cor/t_Z 

    barkas = 1.0 #barkas補正
    LS = 0.0     #LS補正
    delta = 0    #密度補正

    energy_loss  = ((f2)*barkas + LS - delta/2.)*f1 #補正なしのEnergy loss MeV/(g/cm2)
    energy_straggling = f1 * electron_mass * beta2 * gamma * gamma * LSX #Energy straggling  MeV^2/(g/cm2)
    
    return energy_loss, energy_straggling

dedxの定数 dedx_constant=Kを次式のように定義すると、

\displaystyle K=4\pi N_A r_e^2 m_e c^2

諸々の補正項を無視するとdE/dxは次式のように表現される。

\displaystyle \frac{dE}{dx} = K Z_p^2 \frac{Z_t}{A_t}\frac{1}{\beta^2} \left(\ln \left(\frac{2m_e \beta^2 \gamma^2}{I} \right) - \beta^2\right)

energy loss straggling Ωは次式のように表現される。

\displaystyle \frac{d\Omega^2}{dx} = K Z_p^2 \frac{Z_t}{A_t}\frac{1}{\beta^2} m_e \beta^2 \gamma^2 X

345MeV/uの238U、320MeV/uの4Heを、ベリリウムの薄膜に入射したした場合の計算。

print(energy_loss_straggling(p_A=238, p_Z=92, p_Q=92, p_T=345, t_Z=4))
print(energy_loss_straggling(p_A=4, p_Z=2, p_Q=2, p_T=320, t_Z=4))

(22419.411135356077, 1107.0218819306313)
(11.02602078859991, 0.5028755625661142)

と返ってくるはずである。

標的をBe1mmとすると、物質量0.185 g/cm2をかければ、Energy loss と stragglingの2乗の値が得られる。345MeV/uの238U92+ビームのEnergy lossは 4147.6 MeV、Energy stragglingは14.3 MeVとなる。それぞれAで割ると、Energy lossは17.43±0.06 MeV/uとなる。LISE++で出力される値は17.85±0.08 MeV/uである。

320MeV の 4HeビームがBe標的1mmに入射したとき、Energy lossは2.04 MeV、Energy stragglingは0.31 MeVである。それぞれAで割ると、Energy lossは0.510±0.076 MeV/uとなる。LISE++の値は0.508±0.068 MeV/uである。

精度よくEnergy lossを求めるためには、barkas項、delta項、LS項を与える必要がある。

stragglingにかけるべきXの値はATIMAで計算でき、それを与えるとLISE++の値とほぼ一致する。

よくあるStopping powerの図を描く。粒子は陽子、標的は銅(Z=29)にして、横軸はβγにした。

import matplotlib.pyplot as plt
import numpy as np
elosses = []
betagammas = []
for log_energy in range(-10,50):
    energy = 10**(log_energy*0.1)
    eloss, strag = energy_loss_straggling(p_A=1, p_Z=1, p_Q=1, p_T=energy, t_Z=29)
    atomic_mass_unit = 931.4940954 # // 統一原子質量単位 MeV/c^2
    gamma=1.0 + energy/atomic_mass_unit
    beta2=1.0-1.0/(gamma*gamma)
    beta = np.sqrt(beta2)
    betagammas.append(beta*gamma)
    elosses.append(eloss)
plt.rcParams["font.size"] = 14
plt.rcParams['figure.figsize'] = [7, 5]
plt.plot(betagammas,elosses)
plt.xscale("log")
plt.yscale("log")
plt.ylim(1,None)
plt.xlabel("βγ")
plt.ylabel("MeV / (g/cm2)")
plt.show()

PDGの2018年版ブックレットの図33.1と比較する。

概ね問題なさそうだ。

密度補正deltaの部分は比較的簡単にPythonで書けた

def bethek_density_effect(beta, t_Z, beta, gamma):
    '''
    密度効果補正
    '''
    del_0= [
    0.0,  0.00,  0.14, 0.14, 0.14, 0.12, 0.00, 0.00, 0.00, 0.00,
    0.08, 0.08,  0.12, 0.14, 0.14, 0.14, 0.00, 0.00, 0.10, 0.14,
    0.10, 0.12,  0.14, 0.14, 0.14, 0.12, 0.12, 0.10, 0.08, 0.08,
    0.14, 0.14,  0.08, 0.10, 0.00, 0.00, 0.14, 0.14, 0.14, 0.14,
    0.14, 0.14,  0.14, 0.14, 0.14, 0.14, 0.14, 0.14, 0.14, 0.14,
    0.14, 0.14,  0.00, 0.00, 0.14, 0.14, 0.14, 0.14, 0.14, 0.14,
    0.14, 0.14,  0.14, 0.14, 0.14, 0.14, 0.14, 0.14, 0.14, 0.14,
    0.14, 0.14,  0.14, 0.14, 0.08, 0.10, 0.10, 0.12, 0.14, 0.14,
    0.14, 0.14,  0.14, 0.14, 0.00, 0.00, 0.00, 0.14, 0.14, 0.14,
    0.14, 0.14,  0.14, 0.14, 0.14, 0.14, 0.14 ]
    x0= [
    1.8639, 2.2017,  0.1304, 0.0592, 0.0305, -.0178, 1.7378, 1.7541, 1.8433, 2.0735,
    0.2880, 0.1499,  0.1708, 0.2014, 0.1696, 0.1580, 1.5555, 1.7635, 0.3851, 0.3228,
    0.1640, 0.0957,  0.0691, 0.0340, 0.0447, -.0012, -.0187, -.0566, -.0254, 0.0049,
    0.2267, 0.3376,  0.1767, 0.2258, 1.5262, 1.7158, 0.5737, 0.4585, 0.3608, 0.2957,
    0.1785, 0.2267,  0.0949, 0.0599, 0.0576, 0.0563, 0.0657, 0.1281, 0.2406, 0.2879,
    0.3189, 0.3296,  0.0549, 1.5630, 0.5473, 0.4190, 0.3161, 0.2713, 0.2333, 0.1984,
    0.1627, 0.1520,  0.1888, 0.1058, 0.0947, 0.0822, 0.0761, 0.0648, 0.0812, 0.1199,
    0.1560, 0.1965,  0.2117, 0.2167, 0.0559, 0.0891, 0.0819, 0.1484, 0.2021, 0.2756,
    0.3491, 0.3776,  0.4152, 0.4267, 0.4300, 1.5368, 0.6000, 0.5991, 0.4559, 0.4202,
    0.3144, 0.2260,  0.1869, 0.1557, 0.2274, 0.2484, 0.2378]
    gamma = 1/np.sqrt(1-(beta*beta))
    x = np.log(beta * gamma) / 2.3025851
    delta = del_0[t_Z-1] * pow(10.0,(2.*(x-x0[t_Z-1])))
    return delta

全てのコードのまとめ

gitlab.com

参照: C++のCATIMAコード

Python: 4桁の原子量

Python 数学/物理 
standard_weights = [1.008, 4.003,\
6.941, 9.012, 10.81, 12.01, 14.01, 16.00, 19.00, 20.18,\
22.99, 24.31, 26.98, 28.09, 30.97, 32.07, 35.45, 39.95,\
39.10, 40.08, 44.96, 47.87, 50.94, 52.00, 54.94, 55.85, 58.93, 58.69, 63.55, 65.38, 69.72, 72.63, 74.92, 78.97, 79.90, 83.80,\
85.47, 87.62, 88.91, 91.22, 92.91, 95.95, 99, 101.1, 102.9, 106.4, 107.9, 112.4, 114.8, 118.7, 121.8, 127.6, 126.9, 131.3,\
132.9, 137.3,\
138.9, 140.1, 140.9, 144.2, 145, 150.4, 152.0, 157.3, 158.9, 162.5, 164.9, 167.3, 168.9, 173.0, 175.0,\
178.5, 180.9, 183.8, 186.2, 190.2, 192.2, 195.1, 197.0, 200.6, 204.4, 207.2, 209.0, 210, 210, 222,\
223, 226,\
227, 232.0, 231.0, 238.0, 237, 239, 243, 247, 247, 252, 252, 257, 258, 259, 262,\
267, 268, 271, 272, 277, 276, 281, 280, 285, 278, 289, 289, 293, 293, 294]

z2weight = {}
for z,standard_weight in zip(range(1,len(standard_weights)+1),standard_weights):
    z2weight[z]=standard_weight

値は、日本化学会原子量専門委員会が作成した4桁の原子量表(2021)を使った。

こちらもどうぞ

phst.hateblo.jp