qudarap/tests/lhcb_rich_cerenkov_angles_pid.py

0001 #!/usr/bin/env python
0002 """
0003 lhcb_rich_cerenkov_angles_pid.py
0004 ===================================
0005
0006
0007 https://physics.stackexchange.com/questions/471722/cherenkov-light-and-refractive-index
0008
0009
0010 """
0011 import numpy as np, os
0012 import matplotlib.pyplot as plt
0013
0014 def beta(E, mec2):
0015     gamma = E / mec2
0016     return np.sqrt(1 - 1 / np.power(gamma, 2))
0017
0018 def theta(beta, n):
0019     return np.arccos(1 / (beta * n))
0020
0021 mec2_kaon = 493.68  # MeV / c^2
0022 mec2_pion = 139.57  # MeV / c^2
0023 mec2_muon = 105.65  # MeV / c^2
0024
0025 E = np.logspace(np.log10(5e3), 6, 100) # MeV
0026
0027 beta_pion = beta(E, mec2_pion)
0028 beta_kaon = beta(E, mec2_kaon)
0029 beta_muon = beta(E, mec2_muon)
0030
0031 n_1 = 1.0014
0032 n_2 = 1.0005
0033
0034 fig, ax = plt.subplots()
0035
0036 plt.semilogx(E, theta(beta_kaon, n_1), color="crimson", ls="-", label=r"$K, n_1=1.0014$")
0037 plt.semilogx(E, theta(beta_pion, n_1), color="crimson", ls="--", label=r"$\pi, n_1=1.0014$")
0038 plt.semilogx(E, theta(beta_muon, n_1), color="crimson", ls="-.", label=r"$\mu, n_1=1.0014$")
0039
0040 plt.semilogx(E, theta(beta_kaon, n_2), color="k", ls="-", label=r"$K, n_2=1.0005$")
0041 plt.semilogx(E, theta(beta_pion, n_2), color="k", ls="--", label=r"$\pi, n_2=1.0005$")
0042 plt.semilogx(E, theta(beta_muon, n_2), color="k", ls="-.", label=r"$\mu, n_2=1.0005$")
0043
0044 plt.axvline(8e4, color="crimson", ls=":")
0045 plt.axvline(1.1e5, color="k", ls=":")
0046
0047 plt.xlabel("E / MeV")
0048 plt.ylabel(r"$\theta_{\rm C}\,/\,rad$")
0049 plt.legend()
0050 plt.show()
0051
0052 path = "/tmp/qudarap/lhcb_rich_cerenkov_angles_pid.png"
0053 fold = os.path.dirname(path)
0054 if not os.path.isdir(fold):
0055     os.mkdirs(fold)
0056 pass
0057 fig.savefig(path)
0058