Compare commits
4 commits
233684286d
...
d812e60632
| Author | SHA1 | Date | |
|---|---|---|---|
|
d812e60632 | ||
|
|
18e81c76e2 | ||
|
|
80ad5f1693 | ||
|
|
0c6a68b44c |
16 changed files with 1567 additions and 613 deletions
183
Frank-Tamm.py
183
Frank-Tamm.py
|
|
@ -14,9 +14,18 @@ import scienceplots
|
|||
|
||||
from FT_functions import *
|
||||
|
||||
import sys
|
||||
from pathlib import Path
|
||||
|
||||
sys.path.append(str(Path(__file__).resolve().parents[1]))
|
||||
import plotstyle
|
||||
|
||||
MODE = "paper"
|
||||
MODE = "slides"
|
||||
plotstyle.set_style(MODE)
|
||||
|
||||
###Frank-Tamm plot
|
||||
plt.figure(figsize=(15, 10))
|
||||
plt.style.use(['science'])
|
||||
plt.figure()
|
||||
|
||||
m_e = 0.511 # in MeV/c^2
|
||||
m_p = 938.27 # in MeV/c^2
|
||||
|
|
@ -82,8 +91,8 @@ plt.plot(E0*J_to_MeV,n_n, color = 'green', lw=2, label = 'helium')
|
|||
plt.xscale('log')
|
||||
plt.yscale('log')
|
||||
|
||||
plt.ylabel(r'$\frac{dN_{\text{photons}}}{dx} \, \text{in cm}^{-1}$', fontsize=fs+2)
|
||||
plt.xlabel(r'$E_{\text{kin}} \, \text{in MeV}$', fontsize=fs+2)
|
||||
plt.ylabel(r'$\frac{dN_{\text{photons}}}{dx} \, \text{in cm}^{-1}$')
|
||||
plt.xlabel(r'$E_{\text{kin}} \, \text{in MeV}$')
|
||||
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
|
|
@ -96,24 +105,22 @@ plt.axhline(y=57, color='black', linestyle='--', linewidth=2)
|
|||
plt.axhline(y=228, color='green', linestyle='--', linewidth=2)
|
||||
|
||||
plt.xlim(0.5,5000000)
|
||||
plt.ylim(0.7,400)
|
||||
|
||||
plt.title("Number of produced Cherenkov photons in aerogel (300-500nm)",fontsize=fs+3,pad=15)
|
||||
plt.title("Number of produced Cherenkov photons in aerogel (300-500nm)")
|
||||
|
||||
|
||||
#plt.ylim(1,100)
|
||||
plt.legend(fontsize=fs-1, loc="lower right",ncol=1,
|
||||
frameon=True, framealpha=0.4, # leicht transparent
|
||||
fancybox=True,
|
||||
#edgecolor='black', # Rahmenfarbe
|
||||
facecolor='white' # Hintergrundfarbe)
|
||||
)
|
||||
plt.savefig("plots/frank-tamm")
|
||||
#plt.show()
|
||||
# plt.legend(fontsize=fs-1, loc="lower right",ncol=1,
|
||||
# frameon=True, framealpha=0.4, # leicht transparent
|
||||
# fancybox=True,
|
||||
# #edgecolor='black', # Rahmenfarbe
|
||||
# facecolor='white' # Hintergrundfarbe)
|
||||
# )
|
||||
plt.legend(loc="lower right")
|
||||
plotstyle.savefig("frank-tamm", category="Cherenkov")
|
||||
|
||||
|
||||
##################### ANGLE #########################
|
||||
plt.figure(figsize=(15, 10))
|
||||
plt.style.use(['science'])
|
||||
plt.figure()
|
||||
|
||||
# Cherenkov angles
|
||||
theta_e = cherenkov_theta_deg(1.05, E0, m_e)
|
||||
|
|
@ -149,22 +156,18 @@ plt.text(1, theta_max_105, rf'n=1.05: ${theta_max_105:.1f}^\circ$', color='black
|
|||
plt.text(1, theta_max_107, rf'n=1.07: ${theta_max_107:.1f}^\circ$', color='black', fontsize=fs-4, verticalalignment='bottom', horizontalalignment='left')
|
||||
|
||||
plt.xscale('log')
|
||||
plt.xlabel(r'$E_{\mathrm{kin}}\;\text{in MeV}$', fontsize=fs+2)
|
||||
plt.ylabel(r'Cherenkov angle $\theta_C$ in deg', fontsize=fs+2)
|
||||
plt.xlabel(r'$E_{\mathrm{kin}}\;\text{in MeV}$')
|
||||
plt.ylabel(r'Cherenkov angle $\theta_C$ in deg')
|
||||
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
plt.ylim(0,23)
|
||||
plt.tick_params(axis='both', size=7, width=1.5)
|
||||
|
||||
plt.xlim(0.5, 5e6)
|
||||
|
||||
plt.title('Cherenkov angle in aerogel', fontsize=fs+3, pad=15)
|
||||
plt.title('Cherenkov angle in aerogel')
|
||||
|
||||
plt.legend(fontsize=fs-2, loc='lower right',
|
||||
frameon=True, framealpha=0.4, fancybox=True)
|
||||
|
||||
plt.savefig("plots/cherenkov_angle_deg")
|
||||
plt.legend(loc='lower right')
|
||||
plotstyle.savefig("cherenkov_angle_deg", category="Cherenkov")
|
||||
|
||||
|
||||
###############################################################
|
||||
|
|
@ -200,27 +203,18 @@ for l1, l2 in zip(Nlam[:-1], Nlam[1:]):
|
|||
N_lambda = np.array(N_lambda)
|
||||
|
||||
# ---------- Plot ----------
|
||||
plt.figure(figsize=(15, 10))
|
||||
plt.style.use(['science'])
|
||||
plt.figure()
|
||||
|
||||
plt.step(lam_centers, N_lambda, color = 'blue', where='mid',lw=2)
|
||||
#plt.scatter(lam_centers, N_lambda, s=40, color = 'blue')
|
||||
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
plt.grid(alpha=0.3)
|
||||
|
||||
plt.title(
|
||||
r'Cherenkov-Spektrum in aerogel for an $100\,\mathrm{MeV}$ electron',
|
||||
fontsize=fs+3
|
||||
)
|
||||
plt.ylabel(r'$\frac{\mathrm{d}N_{\text{photons}}}{\mathrm{d}x}\,\mathrm{cm}^{-1}$',
|
||||
fontsize=fs + 2)
|
||||
plt.xlabel(r'Wavelength $\lambda$ in nm', fontsize=fs+2)
|
||||
r'Cherenkov-Spektrum in aerogel for an $100\,\mathrm{MeV}$ electron')
|
||||
plt.ylabel(r'$\frac{\mathrm{d}N_{\text{photons}}}{\mathrm{d}x}\,\mathrm{cm}^{-1}$')
|
||||
plt.xlabel(r'Wavelength $\lambda$ in nm')
|
||||
|
||||
plt.xlim(200, 800)
|
||||
plt.tight_layout()
|
||||
plt.savefig("plots/spectrum.pdf")
|
||||
plotstyle.savefig("spectrum", category="Cherenkov")
|
||||
|
||||
|
||||
################## TRANSMITTENZ ###############
|
||||
|
|
@ -250,8 +244,7 @@ T = T_percent / 100.0
|
|||
T_percent_centers = 0.5 * (T_percent[:-1] + T_percent[1:])
|
||||
T_centers = T_percent_centers / 100.0
|
||||
|
||||
fig, ax1 = plt.subplots(figsize=(15, 10))
|
||||
plt.style.use(['science'])
|
||||
fig, ax1 = plt.subplots()
|
||||
|
||||
# --- linke y-Achse: Cherenkov-Spektrum ---
|
||||
ax1.step(
|
||||
|
|
@ -269,10 +262,9 @@ ax1.set_ylabel(
|
|||
fontsize=fs+2, color='blue'
|
||||
)
|
||||
ax1.tick_params(axis='y', labelcolor='blue')
|
||||
ax1.tick_params(axis='both', labelsize=fs)
|
||||
ax1.grid(alpha=0.3)
|
||||
#ax1.tick_params(axis='both', labelsize=fs)
|
||||
ax1.tick_params(axis='y', labelcolor='blue', labelsize=fs)
|
||||
ax1.grid(alpha=0.3)
|
||||
|
||||
|
||||
# --- rechte y-Achse: Transmittanz ---
|
||||
ax2 = ax1.twinx()
|
||||
|
|
@ -294,14 +286,12 @@ ax1.plot(
|
|||
label=r'$\left(\mathrm{d}N/\mathrm{d}x\right)\times T$'
|
||||
)
|
||||
|
||||
ax2.set_ylabel(r'Transmittance (\%)', fontsize=fs+2, color='red')
|
||||
ax2.tick_params(axis='y', labelcolor='red', labelsize=fs)
|
||||
ax2.set_ylabel(r'Transmittance (\%)', color='red')
|
||||
ax2.tick_params(axis='y', labelcolor='red')
|
||||
|
||||
# --- Titel & Limits ---
|
||||
ax1.set_title(
|
||||
r'Effective Cherenkov photon yield in aerogel for a $100\,\mathrm{MeV}$ electron',
|
||||
fontsize=fs+5
|
||||
)
|
||||
r'Effective Cherenkov photon yield in aerogel for a $100\,\mathrm{MeV}$ electron')
|
||||
|
||||
ax1.set_xlim(260, 700)
|
||||
ax1.set_ylim(0, 6.5)
|
||||
|
|
@ -325,24 +315,12 @@ labels = [l_dndx[0], l_T[0], l_dndx[1]]
|
|||
ax1.legend(
|
||||
handles,
|
||||
labels,
|
||||
fontsize=fs,
|
||||
loc='center right',
|
||||
frameon=True,
|
||||
fancybox=True,
|
||||
#edgecolor='black', # Rahmenfarbe
|
||||
facecolor='white' # Hintergrundfarbe
|
||||
loc='center right'
|
||||
)
|
||||
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
|
||||
plt.grid(True, which='major', alpha=0.9, linewidth=1.0)
|
||||
plt.grid(True, which='minor', alpha=0.8, linewidth=0.5)
|
||||
|
||||
#plt.minorticks_on()
|
||||
|
||||
plt.savefig("plots/Transspec.pdf")
|
||||
|
||||
plotstyle.savefig("Transspec", category="Cherenkov")
|
||||
|
||||
################# QUANTEN EFFIZIENZ #####################
|
||||
QE_percent = np.array([
|
||||
|
|
@ -361,33 +339,10 @@ QE_percent = np.array([
|
|||
|
||||
mask = QE_percent[:-1] >= 0
|
||||
|
||||
plt.figure(figsize=(8, 16))
|
||||
plt.style.use(['science'])
|
||||
|
||||
#plt.step(Nlam[:-1][mask],QE_percent[:-1][mask],where='mid',lw=2)
|
||||
plt.plot(lam_centers[mask], QE_percent[:-1][mask], 'x')
|
||||
|
||||
plt.xlabel("Wavelength $\lambda$ in nm", fontsize=fs+2)
|
||||
plt.ylabel("Quantum Efficiency in \%", fontsize=fs+2)
|
||||
plt.title("Photocathode Quantum Efficiency", fontsize=fs+5)
|
||||
plt.xlim(200, 800)
|
||||
plt.yscale('log')
|
||||
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
|
||||
plt.grid(True, which='major', alpha=0.9, linewidth=1.0)
|
||||
plt.grid(True, which='minor', alpha=0.8, linewidth=0.5)
|
||||
|
||||
plt.minorticks_on()
|
||||
plt.savefig("plots/QE.pdf")
|
||||
|
||||
|
||||
####### PREDICTION ##########
|
||||
|
||||
|
||||
fig, ax1 = plt.subplots(figsize=(15, 10))
|
||||
plt.style.use(['science'])
|
||||
fig, ax1 = plt.subplots()
|
||||
|
||||
# --- linke y-Achse: Cherenkov-Spektrum ---
|
||||
ax1.plot(
|
||||
|
|
@ -399,16 +354,13 @@ ax1.plot(
|
|||
)
|
||||
|
||||
|
||||
ax1.set_xlabel(r'Wavelength $\lambda$ in nm', fontsize=fs+2)
|
||||
ax1.set_xlabel(r'Wavelength $\lambda$ in nm')
|
||||
ax1.set_ylabel(
|
||||
r'$\frac{\mathrm{d}N_{\text{photons}}}{\mathrm{d}x}\,\mathrm{cm}^{-1}$',
|
||||
fontsize=fs+2, color='blue'
|
||||
r'$\frac{\mathrm{d}N_{\text{photons}}}{\mathrm{d}x}\,\mathrm{cm}^{-1}$', color='blue'
|
||||
)
|
||||
ax1.tick_params(axis='y', labelcolor='blue')
|
||||
ax1.tick_params(axis='both', labelsize=fs)
|
||||
ax1.grid(alpha=0.3)
|
||||
ax1.tick_params(axis='y', labelcolor='blue', labelsize=fs)
|
||||
ax1.grid(alpha=0.3)
|
||||
ax1.tick_params(axis='both')
|
||||
ax1.tick_params(axis='y', labelcolor='blue')
|
||||
|
||||
# --- rechte y-Achse: Transmittanz ---
|
||||
ax2 = ax1.twinx()
|
||||
|
|
@ -453,14 +405,12 @@ ax1.fill_between(
|
|||
)
|
||||
|
||||
|
||||
ax2.set_ylabel(r'Transmittance (\%)', fontsize=fs+2, color='red')
|
||||
ax2.tick_params(axis='y', labelcolor='red', labelsize=fs)
|
||||
ax2.set_ylabel(r'Transmittance (\%)', color='red')
|
||||
ax2.tick_params(axis='y', labelcolor='red')
|
||||
|
||||
# --- Titel & Limits ---
|
||||
ax1.set_title(
|
||||
r'Effective Cherenkov photon yield in aerogel for a $100\,\mathrm{MeV}$ electron',
|
||||
fontsize=fs+5
|
||||
)
|
||||
r'Effective Cherenkov photon yield in aerogel for a $100\,\mathrm{MeV}$ electron')
|
||||
|
||||
ax1.set_xlim(260, 700)
|
||||
#ax1.set_ylim(0, 6.5)
|
||||
|
|
@ -475,15 +425,30 @@ labels = [l_dndx[0], l_T[0], l_T[1], l_dndx[1], l_dndx[2], l_dndx[3]]
|
|||
ax1.legend(
|
||||
handles,
|
||||
labels,
|
||||
fontsize=fs,
|
||||
loc='upper right',
|
||||
frameon=True,
|
||||
fancybox=True,
|
||||
#framealpha=0.9, # leicht transparent
|
||||
edgecolor='black', # Rahmenfarbe
|
||||
facecolor='white' # Hintergrundfarbe
|
||||
)
|
||||
|
||||
plt.minorticks_on()
|
||||
ax2.set_yscale('log')
|
||||
#ax1.set_yscale('log')
|
||||
|
||||
plotstyle.savefig("prediction_wrong", category="Cherenkov")
|
||||
|
||||
|
||||
####### QE ################
|
||||
|
||||
plt.figure(figsize=(8, 16))
|
||||
|
||||
#plt.step(Nlam[:-1][mask],QE_percent[:-1][mask],where='mid',lw=2)
|
||||
plt.plot(lam_centers[mask], QE_percent[:-1][mask], 'x')
|
||||
|
||||
plt.xlabel("Wavelength $\lambda$ in nm", fontsize=fs+2)
|
||||
plt.ylabel("Quantum Efficiency in \%", fontsize=fs+2)
|
||||
plt.title("Photocathode Quantum Efficiency", fontsize=fs+5)
|
||||
plt.xlim(200, 800)
|
||||
plt.yscale('log')
|
||||
plt.ylim(0.008,50)
|
||||
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
|
||||
|
|
@ -491,10 +456,8 @@ plt.grid(True, which='major', alpha=0.9, linewidth=1.0)
|
|||
plt.grid(True, which='minor', alpha=0.8, linewidth=0.5)
|
||||
|
||||
plt.minorticks_on()
|
||||
ax2.set_yscale('log')
|
||||
#ax1.set_yscale('log')
|
||||
|
||||
#plt.savefig("plots/prediction.pdf")
|
||||
#plt.savefig("plots/QE.pdf")
|
||||
plotstyle.savefig("QE", category="Cherenkov")
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
269
n_dependence.py
269
n_dependence.py
|
|
@ -6,7 +6,15 @@ import scienceplots
|
|||
|
||||
from FT_functions import *
|
||||
|
||||
fs=40
|
||||
import sys
|
||||
from pathlib import Path
|
||||
|
||||
sys.path.append(str(Path(__file__).resolve().parents[1]))
|
||||
import plotstyle
|
||||
|
||||
MODE = "paper"
|
||||
#MODE = "slides"
|
||||
plotstyle.set_style(MODE)
|
||||
|
||||
#### Dependence on n #########
|
||||
|
||||
|
|
@ -34,6 +42,10 @@ z = z_e
|
|||
lam1 = 300
|
||||
lam2 = 500
|
||||
|
||||
# =================================================
|
||||
# Photonenzahl vs n
|
||||
# =================================================
|
||||
|
||||
# Brechungsindex-Bereich
|
||||
n_values = np.linspace(1.04, 1.10, 100)
|
||||
|
||||
|
|
@ -52,34 +64,18 @@ for n_val in n_values:
|
|||
|
||||
N_total = np.array(N_total)
|
||||
|
||||
# ---------- Plot ----------
|
||||
plt.figure(figsize=(15, 10))
|
||||
plt.style.use(['science'])
|
||||
|
||||
plt.plot(n_values, N_total, lw=2, color='blue')
|
||||
|
||||
plt.xlabel(r'Refractive index $n$', fontsize=fs+2)
|
||||
plt.ylabel(
|
||||
r'$\frac{dN_{\mathrm{photons}}}{dx}$ in $\mathrm{cm}^{-1}$',
|
||||
fontsize=fs+2
|
||||
)
|
||||
|
||||
plt.title(
|
||||
r'Cherenkov photon yield'
|
||||
'\n$100\,\mathrm{MeV}$ electron, 300–500 nm',
|
||||
fontsize=fs+3
|
||||
)
|
||||
|
||||
plt.grid(alpha=0.3)
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
|
||||
plt.tight_layout()
|
||||
plt.savefig("plots/photons_n.pdf")
|
||||
plt.figure()
|
||||
plt.plot(n_values, N_total, lw=2, color='blue', label='electron')
|
||||
plt.xlabel(r'Refractive index $n$')
|
||||
plt.ylabel(r'$\frac{dN_{\mathrm{photons}}}{dx}$ in $\mathrm{cm}^{-1}$')
|
||||
plt.title(r'Cherenkov photon yield ($100\,\mathrm{MeV}$ electron)')
|
||||
plt.legend()
|
||||
plotstyle.savefig("photons_n", category="Cherenkov")
|
||||
|
||||
|
||||
######### PHOTON EMISSSION PRO WELLENLÄNGE ###############################
|
||||
|
||||
# =================================================
|
||||
# Spektrum vs n
|
||||
# =================================================
|
||||
|
||||
# ---------- Wellenlängenbins ----------
|
||||
Nlam = np.arange(200, 810, 10) # nm
|
||||
|
|
@ -95,8 +91,7 @@ colors = {
|
|||
1.09: 'tab:red'
|
||||
}
|
||||
|
||||
plt.figure(figsize=(15, 10))
|
||||
plt.style.use(['science'])
|
||||
plt.figure()
|
||||
|
||||
for n_val in n_list:
|
||||
N_lambda = []
|
||||
|
|
@ -122,29 +117,17 @@ for n_val in n_list:
|
|||
label=fr'$n={n_val:.2f}$'
|
||||
)
|
||||
|
||||
# ---------- Plot-Format ----------
|
||||
plt.xlabel(r'Wavelength $\lambda$ in nm', fontsize=fs+2)
|
||||
plt.ylabel(
|
||||
r'$\frac{dN_{\mathrm{photons}}}{dx}$ in $\mathrm{cm}^{-1}$',
|
||||
fontsize=fs+2
|
||||
)
|
||||
|
||||
plt.title(
|
||||
'Cherenkov spectrum in aerogel ($100\,\mathrm{MeV}$ electron)',
|
||||
fontsize=fs+3
|
||||
)
|
||||
|
||||
plt.xlabel(r'Wavelength $\lambda$ in nm')
|
||||
plt.ylabel(r'$\frac{dN_{\mathrm{photons}}}{dx}$ in $\mathrm{cm}^{-1}$')
|
||||
plt.title(r'Cherenkov spectrum in aerogel ($100\,\mathrm{MeV}$ electron)')
|
||||
plt.xlim(200, 800)
|
||||
plt.grid(alpha=0.3)
|
||||
plt.legend(fontsize=fs-2, frameon=True)
|
||||
plt.legend()
|
||||
plotstyle.savefig("spectrum_n", category="Cherenkov")
|
||||
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
# =================================================
|
||||
# Spektrum × Transmittanz
|
||||
# =================================================
|
||||
|
||||
plt.tight_layout()
|
||||
plt.savefig("plots/spectrum_n.pdf")
|
||||
|
||||
################# SPEKTRUM UND TRANSMITTENZ ###########
|
||||
T_percent = np.array([
|
||||
0, 0, 0, 0.2, 1.1, 2.5, 4.5, # 200–260
|
||||
7.1, 10.4, 14.4, 18.8, # 270–300
|
||||
|
|
@ -171,9 +154,7 @@ T = T_percent / 100.0
|
|||
T_percent_centers = 0.5 * (T_percent[:-1] + T_percent[1:])
|
||||
T_centers = T_percent_centers / 100.0
|
||||
|
||||
fig, ax1 = plt.subplots(figsize=(15, 10))
|
||||
plt.style.use(['science'])
|
||||
|
||||
fig, ax1 = plt.subplots()
|
||||
# ---------- Cherenkov-Spektren ----------
|
||||
spectra = {}
|
||||
|
||||
|
|
@ -211,62 +192,23 @@ for n_val in n_list:
|
|||
#label=fr'$n={n_val:.2f},\ (\mathrm{{d}}N/\mathrm{{d}}x)\times T$'
|
||||
)
|
||||
|
||||
# ---------- Linke Achse ----------
|
||||
ax1.set_xlabel(r'Wavelength $\lambda$ in nm', fontsize=fs+2)
|
||||
ax1.set_ylabel(
|
||||
r'$\frac{\mathrm{d}N_{\mathrm{photons}}}{\mathrm{d}x}$ in $\mathrm{cm}^{-1}$',
|
||||
fontsize=fs+2
|
||||
)
|
||||
ax1.tick_params(axis='both', labelsize=fs)
|
||||
ax1.grid(True, which='major', alpha=0.8, linewidth=1.0)
|
||||
ax1.grid(True, which='minor', alpha=0.5, linewidth=0.5)
|
||||
ax1.minorticks_on()
|
||||
|
||||
ax1.set_xlabel(r'Wavelength $\lambda$ in nm')
|
||||
ax1.set_ylabel(r'$\frac{dN_{\mathrm{photons}}}{dx}$ in $\mathrm{cm}^{-1}$')
|
||||
ax1.set_xlim(260, 700)
|
||||
ax1.set_ylim(0, 11)
|
||||
# ---------- Rechte Achse: Transmittanz ----------
|
||||
|
||||
ax2 = ax1.twinx()
|
||||
ax2.plot(lam_centers[mask], T_percent[:-1][mask], lw=2, linestyle='--', color='black', label='Transmittance')
|
||||
ax2.set_ylabel(r'Transmittance in (\%)')
|
||||
|
||||
ax2.plot(
|
||||
lam_centers[mask],
|
||||
T_percent[:-1][mask],
|
||||
linestyle='--',
|
||||
color='black',
|
||||
lw=2,
|
||||
label='Transmittance'
|
||||
)
|
||||
|
||||
ax2.set_ylabel(r'Transmittance in (\%)', fontsize=fs+2)
|
||||
ax2.tick_params(axis='y', labelsize=fs)
|
||||
|
||||
# ---------- Titel ----------
|
||||
ax1.set_title(
|
||||
r'Cherenkov spectrum and effective photon yield in aerogel'
|
||||
'\n$100\,\mathrm{MeV}$ electron',
|
||||
fontsize=fs+4
|
||||
)
|
||||
|
||||
# ---------- Gemeinsame Legende ----------
|
||||
h1, l1 = ax1.get_legend_handles_labels()
|
||||
h2, l2 = ax2.get_legend_handles_labels()
|
||||
ax1.legend(h1 + h2, l1 + l2, loc='center right')
|
||||
plotstyle.savefig("Transspec_multi_n", category="Cherenkov")
|
||||
|
||||
ax1.legend(
|
||||
h1 + h2,
|
||||
l1 + l2,
|
||||
fontsize=fs-2,
|
||||
loc='upper right',
|
||||
frameon=True,
|
||||
fancybox=True,
|
||||
framealpha=0.85,
|
||||
facecolor='white'
|
||||
)
|
||||
|
||||
plt.tight_layout()
|
||||
plt.savefig("plots/Transspec_multi_n.pdf")
|
||||
|
||||
|
||||
|
||||
################### Vergleichbare Phototnen mit Transmittenz ####
|
||||
# =================================================
|
||||
# Effektiver Photonenzahl
|
||||
# =================================================
|
||||
|
||||
lam1 = 300
|
||||
lam2 = 500
|
||||
|
|
@ -321,30 +263,16 @@ for n_val in n_values:
|
|||
|
||||
N_eff_total = np.array(N_eff_total)
|
||||
|
||||
plt.figure(figsize=(15, 10))
|
||||
plt.style.use(['science'])
|
||||
plt.plot( n_values, N_eff_total, lw=2, color='black' )
|
||||
plt.xlabel(r'Refractive index $n$', fontsize=fs+2)
|
||||
plt.ylabel(
|
||||
r'Effective Cherenkov photon yield'
|
||||
'\n$\int_{300}^{500}\! (\mathrm{d}N/\mathrm{d}x)\,T(\lambda)\,\mathrm{d}\lambda'
|
||||
r'\;\text{in}\;\mathrm{cm}^{-1}$',
|
||||
fontsize=fs+2
|
||||
)
|
||||
plt.title( r'Effective Cherenkov photon yield vs. refractive index' '\n$100\,\mathrm{MeV}$ electron, including aerogel transmittance', fontsize=fs+3 )
|
||||
plt.grid(True, which='major', alpha=0.8, linewidth=1.0)
|
||||
plt.grid(True, which='minor', alpha=0.5, linewidth=0.5)
|
||||
plt.minorticks_on()
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
plt.tight_layout()
|
||||
plt.savefig("plots/photons_n_effective.pdf")
|
||||
plt.figure()
|
||||
plt.plot(n_values, N_eff_total, lw=2, color='black')
|
||||
plt.xlabel(r'Refractive index $n$')
|
||||
plt.ylabel(r'Effective Cherenkov photon yield ')
|
||||
plt.title(r'Cherenkov photon yield (300–500 nm) vs. refractive index')
|
||||
plotstyle.savefig("photons_n_effective", category="Cherenkov")
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
##### Schwellenenergien ########
|
||||
# =================================================
|
||||
# Schwellenenergie für alle Teilchen
|
||||
# =================================================
|
||||
|
||||
E_thr_e = threshold_energy(n_values, m_e) * J_to_MeV
|
||||
E_thr_p = threshold_energy(n_values, m_p) * J_to_MeV
|
||||
|
|
@ -352,40 +280,21 @@ E_thr_mu = threshold_energy(n_values, m_mu) * J_to_MeV
|
|||
E_thr_a = threshold_energy(n_values, m_a) * J_to_MeV
|
||||
|
||||
|
||||
# ---------- Plot ----------
|
||||
plt.figure(figsize=(15, 10))
|
||||
plt.style.use(['science'])
|
||||
|
||||
plt.plot(n_values, E_thr_e, lw=2, color='blue', label='electron')
|
||||
plt.plot(n_values, E_thr_p, lw=2, color='red', label='proton')
|
||||
plt.plot(n_values, E_thr_a, lw=2, color='green', label='helium')
|
||||
plt.plot(n_values, E_thr_mu, lw=2, color=(0.5, 0.3, 0.8), label='muon')
|
||||
|
||||
plt.figure()
|
||||
plt.plot(n_values, E_thr_e, color='blue', lw=2, label='electron')
|
||||
plt.plot(n_values, E_thr_p, color='red', lw=2, label='proton')
|
||||
plt.plot(n_values, E_thr_a, color='green', lw=2, label='helium')
|
||||
plt.plot(n_values, E_thr_mu, color=(0.5, 0.3, 0.8), lw=2, label='muon')
|
||||
plt.yscale('log')
|
||||
plt.xlabel(r'Refractive index $n$')
|
||||
plt.ylabel(r'$E_{\mathrm{thr}}$ in MeV')
|
||||
plt.title(r'Cherenkov threshold energy')
|
||||
plt.legend()
|
||||
plotstyle.savefig("threshold_energy_n", category="Cherenkov")
|
||||
|
||||
plt.xlabel(r'Refractive index $n$', fontsize=fs+2)
|
||||
plt.ylabel(r'$E_{\mathrm{thr}}$ in MeV', fontsize=fs+2)
|
||||
|
||||
plt.title(
|
||||
r'Cherenkov threshold energy',
|
||||
fontsize=30,
|
||||
pad=15
|
||||
)
|
||||
|
||||
plt.grid(alpha=0.3)
|
||||
plt.legend(
|
||||
fontsize=fs-2,
|
||||
frameon=True,
|
||||
framealpha=0.4,
|
||||
fancybox=True
|
||||
)
|
||||
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
plt.minorticks_on()
|
||||
plt.tight_layout()
|
||||
plt.savefig("plots/threshold_energy_n.pdf")
|
||||
|
||||
# =================================================
|
||||
# Schwellenenergie Proton (Zoom)
|
||||
# =================================================
|
||||
|
||||
# ---------- n-Bereich ----------
|
||||
n_values_p = np.linspace(1.05, 1.10, 200)
|
||||
|
|
@ -393,47 +302,13 @@ n_values_p = np.linspace(1.05, 1.10, 200)
|
|||
# ---------- Schwellenenergie Proton ----------
|
||||
E_thr_p = threshold_energy(n_values_p, m_p) * J_to_MeV
|
||||
|
||||
# ---------- Plot ----------
|
||||
plt.figure(figsize=(15, 10))
|
||||
plt.style.use(['science'])
|
||||
|
||||
plt.plot(
|
||||
n_values_p,
|
||||
E_thr_p/1000,
|
||||
lw=2,
|
||||
color='red',
|
||||
label='proton'
|
||||
)
|
||||
|
||||
#plt.yscale('log')
|
||||
|
||||
plt.xlabel(r'Refractive index $n$', fontsize=fs+2)
|
||||
plt.ylabel(
|
||||
r'$E_{\mathrm{thr}}$ in GeV',
|
||||
fontsize=fs+2
|
||||
)
|
||||
|
||||
plt.title(
|
||||
r'Cherenkov threshold energy for protons'
|
||||
'\n$1.05 \le n \le 1.10$',
|
||||
fontsize=fs+5,
|
||||
pad=15
|
||||
)
|
||||
|
||||
plt.grid(alpha=0.3)
|
||||
plt.legend(
|
||||
fontsize=fs-2,
|
||||
frameon=True,
|
||||
framealpha=0.4,
|
||||
fancybox=True
|
||||
)
|
||||
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
plt.grid(which='minor', linestyle='--', linewidth=0.5, alpha=0.5)
|
||||
plt.minorticks_on()
|
||||
plt.tight_layout()
|
||||
plt.savefig("plots/threshold_energy_proton.pdf")
|
||||
plt.figure()
|
||||
plt.plot(n_values_p, E_thr_p/1000, color='red', lw=2, label='proton')
|
||||
plt.xlabel(r'Refractive index $n$')
|
||||
plt.ylabel(r'$E_{\mathrm{thr}}$ in GeV')
|
||||
plt.title(r'Cherenkov threshold energy for protons')
|
||||
plt.legend()
|
||||
plotstyle.savefig("threshold_energy_proton", category="Cherenkov")
|
||||
|
||||
|
||||
plt.show()
|
||||
|
|
@ -2,7 +2,16 @@ import matplotlib.pyplot as plt
|
|||
import numpy as np
|
||||
from scipy.integrate import simpson, quad
|
||||
import scienceplots
|
||||
plt.style.use(['science'])
|
||||
import sys
|
||||
from pathlib import Path
|
||||
|
||||
sys.path.append(str(Path(__file__).resolve().parents[1]))
|
||||
import plotstyle
|
||||
|
||||
MODE = "paper"
|
||||
#MODE = "slides"
|
||||
plotstyle.set_style(MODE)
|
||||
|
||||
color_cycle = plt.rcParams['axes.prop_cycle'].by_key()['color']
|
||||
|
||||
def FFS(T=np.logspace(-2,2,100), a=0.4): # Force Field Solution
|
||||
|
|
@ -28,7 +37,7 @@ x2l, x2h, y2 = np.loadtxt('datafiles/adriani-etal-2015-electron-spectra.dat', sk
|
|||
|
||||
#print(np.sqrt(x1l*x1h), x1l, x1h)
|
||||
|
||||
figsize = (8,5)
|
||||
#figsize = (8,5)
|
||||
linewidth = 12
|
||||
# fig, ax = plt.subplots(figsize=(6,3.5))
|
||||
|
||||
|
|
@ -352,7 +361,7 @@ linewidth = 12
|
|||
|
||||
|
||||
|
||||
fig, ax = plt.subplots(figsize=figsize)
|
||||
fig, ax = plt.subplots()
|
||||
|
||||
#energy ranges (xmin and xmax are not in data coordinates!!!)
|
||||
xlow, xhigh = 1e-2, 1e2
|
||||
|
|
@ -360,6 +369,8 @@ def scale(x):
|
|||
return (x-xlow)/(xhigh-xlow)
|
||||
ax.set_xlim(xlow,xhigh)
|
||||
|
||||
linewidth=22
|
||||
|
||||
ax.plot(np.sqrt(x1l*x1h), y1, c='k', ls='none', label='Protons p$^+$', marker = 'o',zorder=3)
|
||||
ax.plot(np.sqrt(x2l*x2h), y2/1000, c='k', ls='none', label='Electrons e$^-$', marker = 's',zorder=3)
|
||||
ax.hlines(y=[4e-1],xmin=10e-3,xmax=107e-3,color='tab:orange',alpha=0.5,linewidth=linewidth,label='Range p$^+$')
|
||||
|
|
@ -378,22 +389,25 @@ ax.hlines(y=[2.5e-4],xmin=1.0,xmax=11e3,color='tab:blue',alpha=0.2,linewidth=lin
|
|||
# ax.plot(T_a/1000, sp*1.e4, c='b', ls='-', label='FFS')
|
||||
# ax.plot(T_a/1000, se*1.e4, c='b', ls='-')
|
||||
|
||||
plt.text(0.91, 0.745, 'HET', fontsize=15, rotation=0, transform=plt.gcf().transFigure, va='center')
|
||||
plt.text(0.91, 0.575, 'AMS', fontsize=15, rotation=0, transform=plt.gcf().transFigure, va='center')
|
||||
plt.text(0.91, 0.41, 'KET', fontsize=15, rotation=0, transform=plt.gcf().transFigure, va='center')
|
||||
plt.text(0.91, 0.24, 'ProHEPaM', fontsize=15, rotation=0, transform=plt.gcf().transFigure, va='center')
|
||||
plt.text(0.91, 0.745, 'HET', fontsize=20, rotation=0, transform=plt.gcf().transFigure, va='center')
|
||||
plt.text(0.91, 0.575, 'AMS', fontsize=20, rotation=0, transform=plt.gcf().transFigure, va='center')
|
||||
plt.text(0.91, 0.41, 'KET', fontsize=20, rotation=0, transform=plt.gcf().transFigure, va='center')
|
||||
plt.text(0.91, 0.24, 'ProHEPaM', fontsize=20, rotation=0, transform=plt.gcf().transFigure, va='center')
|
||||
|
||||
ax.set_ylim(5e-5,2e0)
|
||||
ax.set_xscale('log')
|
||||
ax.set_yscale('log')
|
||||
ax.set_xlabel(r'Kinetic Energy $E_{\rm kin}$ in GeV',fontsize=14)
|
||||
ax.set_ylabel(r'Differential Flux $\Phi$ in (m$^2$ s sr MeV)$^{-1}$', fontsize=14)
|
||||
ax.set_xlabel(r'Kinetic Energy $E_{\rm kin}$ in GeV')
|
||||
ax.set_ylabel(r'Differential Flux $\Phi$ in (m$^2$ s sr MeV)$^{-1}$')
|
||||
# ax.set_xlabel(r'$E_{\rm kin}$ in GeV',fontsize=14)
|
||||
# ax.set_ylabel(r'd$J$/d$E$ in (m$^2$ s sr MeV)$^{-1}$', fontsize=14)
|
||||
|
||||
ax.grid(visible=True, which='both', axis='both', ls='--')
|
||||
plt.legend(fontsize=12,loc='upper right')
|
||||
plt.title('GCR Energy Spectra',size=16)
|
||||
plt.savefig('plots/adriani-etal-combined-e-p_all_prohepam.pdf')
|
||||
plt.legend(loc='upper right')
|
||||
plt.title('GCR Energy Spectra')
|
||||
plotstyle.savefig('adriani-etal-combined-e-p_all_prohepam', category="other")
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
# p_flux = simpson(y1*1000, np.sqrt(x1l*x1h))#, dx = x1h - x1l)
|
||||
|
|
|
|||
|
|
@ -69,10 +69,11 @@ plt.xlabel("Depth z in cm", fontsize=fs)
|
|||
plt.ylabel(r"$\langle \mathrm{d}E/\mathrm{d}z \rangle$ in MeV/cm", fontsize=fs)
|
||||
plt.title("Longitudinal energy deposition in BGO (Geant4)", fontsize=fs+2)
|
||||
plt.xlim(0, z_max_cm)
|
||||
plt.grid(True, ls="--", lw=0.6)
|
||||
plt.grid(True)
|
||||
plt.tick_params(labelsize=fs-2)
|
||||
plt.legend(title="Initial Energy", fontsize=fs-4, title_fontsize=fs-2, frameon=True, framealpha=0.85, facecolor="white")
|
||||
|
||||
plt.tight_layout()
|
||||
plt.savefig("plots/G4_longitudinal_profile.png", dpi=300)
|
||||
plt.show()
|
||||
plt.savefig("plots/G4_longitudinal_profile.pdf")
|
||||
plt.show()
|
||||
|
|
|
|||
202
shower_long_comparison.py
Normal file
202
shower_long_comparison.py
Normal file
|
|
@ -0,0 +1,202 @@
|
|||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Longitudinal shower profile comparison
|
||||
Theory (analytic, dashed) vs Geant4 simulation (steps)
|
||||
With residuals panel
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
from scipy.stats import gamma
|
||||
from matplotlib.lines import Line2D
|
||||
|
||||
from pathlib import Path
|
||||
import sys
|
||||
sys.path.append(str(Path(__file__).resolve().parents[1]))
|
||||
import plotstyle
|
||||
|
||||
import scienceplots # optional, falls noch andere Style-Features gebraucht werden
|
||||
|
||||
# =====================================================
|
||||
# Plot Style
|
||||
# =====================================================
|
||||
MODE = "slides" # "paper" oder "slides"
|
||||
plotstyle.set_style(MODE)
|
||||
|
||||
# -------------------------------------------------
|
||||
# Input files (Geant4)
|
||||
# -------------------------------------------------
|
||||
files = [
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_100MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_200MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_500MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_1GeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_2GeV_0.hits",
|
||||
]
|
||||
|
||||
energies = [100, 200, 500, 1000, 2000] # MeV
|
||||
labels = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
|
||||
colors = ["C0", "C1", "C2", "C3", "C4"]
|
||||
|
||||
# -------------------------------------------------
|
||||
# Material (BGO)
|
||||
# -------------------------------------------------
|
||||
X0 = 1.118 # cm
|
||||
Ec = 10.5 # MeV
|
||||
b = 0.5
|
||||
|
||||
# -------------------------------------------------
|
||||
# Plot parameters
|
||||
# -------------------------------------------------
|
||||
z_max_cm = 10.0
|
||||
n_bins = 100
|
||||
|
||||
# -------------------------------------------------
|
||||
# Theory functions
|
||||
# -------------------------------------------------
|
||||
def tmax(E):
|
||||
return np.log(E / Ec) - 0.5
|
||||
|
||||
def dEdz_theory(z, E):
|
||||
t = z / X0
|
||||
a = 1 + b * tmax(E)
|
||||
return E * gamma.pdf(t, a, scale=1 / b) / X0 # MeV/cm/nuc
|
||||
|
||||
# =====================================================
|
||||
# Figure layout
|
||||
# =====================================================
|
||||
fig, (ax, ax_res) = plt.subplots(
|
||||
2, 1,
|
||||
figsize=(plt.rcParams["figure.figsize"][0], plt.rcParams["figure.figsize"][1]),
|
||||
sharex=True,
|
||||
gridspec_kw={"height_ratios": [3, 1], "hspace": 0.08}
|
||||
)
|
||||
|
||||
energy_handles = []
|
||||
|
||||
# -------------------------------------------------
|
||||
# Loop energies
|
||||
# -------------------------------------------------
|
||||
for file, E, label, color in zip(files, energies, labels, colors):
|
||||
|
||||
# -------------------------
|
||||
# Geant4 simulation
|
||||
# -------------------------
|
||||
df = pd.read_csv(file, sep="\t")
|
||||
N_events = df["event"].nunique()
|
||||
|
||||
z_edges = np.linspace(0.0, z_max_cm, n_bins + 1)
|
||||
dz = z_edges[1] - z_edges[0]
|
||||
z_centers = 0.5 * (z_edges[:-1] + z_edges[1:])
|
||||
|
||||
E_sum, _ = np.histogram(df["z"], bins=z_edges, weights=df["edep"])
|
||||
profile_sim = E_sum / (N_events * dz)
|
||||
|
||||
ax.step(
|
||||
z_centers,
|
||||
profile_sim,
|
||||
where="mid",
|
||||
lw=plt.rcParams["lines.linewidth"],
|
||||
color=color,
|
||||
alpha=0.9
|
||||
)
|
||||
|
||||
# -------------------------
|
||||
# Theory
|
||||
# -------------------------
|
||||
z_th = np.linspace(0, z_max_cm, 600)
|
||||
profile_th = dEdz_theory(z_th, E)
|
||||
|
||||
ln, = ax.plot(
|
||||
z_th,
|
||||
profile_th,
|
||||
color=color,
|
||||
linewidth=plt.rcParams["lines.linewidth"],
|
||||
linestyle="--"
|
||||
)
|
||||
|
||||
energy_handles.append(Line2D([0], [0], color=color, lw=plt.rcParams["lines.linewidth"], linestyle="-"))
|
||||
|
||||
# -------------------------
|
||||
# Residuals
|
||||
# -------------------------
|
||||
th_interp = np.interp(z_centers, z_th, profile_th)
|
||||
mask = th_interp > 1e-6
|
||||
|
||||
residuals = np.zeros_like(profile_sim)
|
||||
residuals[mask] = profile_sim[mask]-th_interp[mask]
|
||||
|
||||
ax_res.step(
|
||||
z_centers,
|
||||
residuals,
|
||||
color=color,
|
||||
linewidth=plt.rcParams["lines.linewidth"]/1.5
|
||||
)
|
||||
|
||||
# schwarze Linie bei y=0
|
||||
ax_res.axhline(0.0, color="black", linestyle=":", linewidth=plt.rcParams["lines.linewidth"]/1.5)
|
||||
|
||||
# -------------------------------------------------
|
||||
# Main axis
|
||||
# -------------------------------------------------
|
||||
ax.set_ylabel(r"$\langle \mathrm{d}E/\mathrm{d}z \rangle$ in MeV/cm", fontsize=plt.rcParams["axes.labelsize"])
|
||||
ax.set_title("Longitudinal Energy Deposition in BGO", fontsize=plt.rcParams["axes.titlesize"], pad=plt.rcParams.get("axes.titlepad", 15))
|
||||
ax.set_xlim(0, z_max_cm)
|
||||
ax.grid(True)
|
||||
ax.tick_params(labelsize=plt.rcParams["xtick.labelsize"])
|
||||
|
||||
# -------------------------------------------------
|
||||
# Residual axis
|
||||
# -------------------------------------------------
|
||||
ax_res.set_xlabel("Depth z in cm", fontsize=plt.rcParams["axes.labelsize"])
|
||||
ax_res.set_ylabel(r"$\mathrm{G4sim} - \mathrm{Theory}$", fontsize=plt.rcParams["ytick.labelsize"])
|
||||
ax_res.axhline(0.0, color="black", linestyle=":", linewidth=plt.rcParams["lines.linewidth"]/1.5)
|
||||
ax_res.grid(True, which="both", linestyle="--", linewidth=0.4)
|
||||
ax_res.tick_params(labelsize=plt.rcParams["xtick.labelsize"])
|
||||
ax_res.set_xlim(0, z_max_cm)
|
||||
ax_res.set_ylim(-20,20)
|
||||
|
||||
# -------------------------------------------------
|
||||
# Legends
|
||||
# -------------------------------------------------
|
||||
# Energy legend (colors only)
|
||||
leg1 = ax.legend(
|
||||
handles=energy_handles,
|
||||
labels=labels,
|
||||
title="Initial Energy",
|
||||
fontsize=plt.rcParams["legend.fontsize"]-2,
|
||||
title_fontsize=plt.rcParams["legend.fontsize"],
|
||||
loc="upper right",
|
||||
frameon=True,
|
||||
framealpha=plt.rcParams.get("legend.framealpha", 0.85),
|
||||
facecolor="white"
|
||||
)
|
||||
|
||||
# Style legend
|
||||
style_handles = [
|
||||
Line2D([0], [0], color="black", linewidth=plt.rcParams["lines.linewidth"], linestyle="--"),
|
||||
Line2D([0], [0], color="black", linewidth=plt.rcParams["lines.linewidth"], linestyle="-"),
|
||||
]
|
||||
|
||||
style_labels = ["Theory", "Geant4 simulation"]
|
||||
|
||||
leg2 = ax.legend(
|
||||
handles=style_handles,
|
||||
labels=style_labels,
|
||||
fontsize=plt.rcParams["legend.fontsize"]-3,
|
||||
loc="upper left",
|
||||
frameon=True,
|
||||
framealpha=plt.rcParams.get("legend.framealpha", 0.85),
|
||||
facecolor="white"
|
||||
)
|
||||
|
||||
ax.add_artist(leg1)
|
||||
|
||||
# -------------------------------------------------
|
||||
# Save
|
||||
# -------------------------------------------------
|
||||
plt.tight_layout()
|
||||
plotstyle.savefig("shower_long_comparison", category="BGO")
|
||||
plt.show()
|
||||
184
shower_map.py
184
shower_map.py
|
|
@ -1,7 +1,8 @@
|
|||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Overview: 2D map, transversal and longitudinal shower profile theory
|
||||
2D map, transversal and longitudinal shower profile theory (BGO, NKG model)
|
||||
Paper-style version, fully cleaned
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
|
|
@ -9,153 +10,150 @@ import matplotlib.pyplot as plt
|
|||
from matplotlib.gridspec import GridSpec
|
||||
from scipy.stats import gamma
|
||||
from matplotlib.lines import Line2D
|
||||
import scienceplots
|
||||
from pathlib import Path
|
||||
import sys
|
||||
import scienceplots
|
||||
|
||||
# Plotstyle importieren
|
||||
sys.path.append(str(Path(__file__).resolve().parents[1]))
|
||||
import plotstyle
|
||||
|
||||
# ============================================================
|
||||
# Paper-style setzen
|
||||
# ============================================================
|
||||
plotstyle.set_style("paper")
|
||||
|
||||
# --------------------------
|
||||
# Funktionen
|
||||
# Physikalische Funktionen
|
||||
# --------------------------
|
||||
def tmax(E, Ec):
|
||||
return np.log(E / Ec) - 0.5
|
||||
|
||||
def longitudinal_profile(t, E, Ec, b=0.5):
|
||||
t_max_val = tmax(E, Ec)
|
||||
a = b * t_max_val + 1
|
||||
a = b * tmax(E, Ec) + 1
|
||||
return gamma.pdf(t, a, scale=1/b) * E
|
||||
|
||||
def dEdz(z, E, X0=1.118, b=0.5, Ec=10.5):
|
||||
t = z / X0
|
||||
a = 1 + b * (np.log(E/Ec) - 0.5)
|
||||
return E * gamma.pdf(t, a, scale=1/b) / X0 # MeV/cm/nuc
|
||||
return E * gamma.pdf(z/X0, a, scale=1/b) / X0
|
||||
|
||||
def rho_r(r, Rm=2.26):
|
||||
return np.exp(-r**2/(2*Rm**2)) / (2*np.pi*Rm**2)
|
||||
def rho_r(r, Rm_cm=2.26, s=1.0):
|
||||
x = r / Rm_cm
|
||||
core = np.zeros_like(r)
|
||||
mask = r > 0
|
||||
core[mask] = (x[mask]**(s-2)) * ((1 + x[mask])**(s-4.5))
|
||||
if r[0] == 0:
|
||||
core[0] = ((1e-9 / Rm_cm)**(s-2)) * ((1 + 1e-9 / Rm_cm)**(s-4.5))
|
||||
norm = np.trapz(2*np.pi*r*core, r)
|
||||
return core / norm
|
||||
|
||||
# --------------------------
|
||||
# Parameter
|
||||
# Parameter (BGO)
|
||||
# --------------------------
|
||||
X0_cm = 1.118
|
||||
Ec_MeV = 10.5
|
||||
b = 0.5
|
||||
Rm_cm = 2.26
|
||||
max_depth_cm = 10
|
||||
b = 0.5
|
||||
|
||||
energies_GeV = [0.1, 0.2, 0.5, 1.0, 2.0]
|
||||
energies_MeV = [E*1000 for E in energies_GeV]
|
||||
colors = ['C0', 'C1', 'C2', 'C3', 'C4']
|
||||
|
||||
fs = 22 # Schriftgröße
|
||||
|
||||
# Heatmap Konturen in Prozent
|
||||
heatmap_levels = [1, 10, 50, 90]
|
||||
linestyles = {1:':', 10:'--', 50:'-.', 90:'-'}
|
||||
heatmap_levels = [0.1, 1, 10]
|
||||
linestyles = [":", "--", "-"]
|
||||
|
||||
# --------------------------
|
||||
# Gitter
|
||||
# --------------------------
|
||||
r = np.linspace(0, 8, 300)
|
||||
r_max = 2.5
|
||||
r = np.linspace(0, r_max, 300)
|
||||
z = np.linspace(0, max_depth_cm, 400)
|
||||
Z, R = np.meshgrid(z, r, indexing="ij")
|
||||
|
||||
# --------------------------
|
||||
# Figure und GridSpec
|
||||
# Figure und Layout
|
||||
# --------------------------
|
||||
fig = plt.figure(figsize=(14,10))
|
||||
plt.style.use(['science'])
|
||||
fig.suptitle("Electromagnetic shower in BGO (theory)", fontsize=fs+2, y=0.97)
|
||||
gs = GridSpec(2,2, width_ratios=[4,1], height_ratios=[4,1], hspace=0.3, wspace=0.1)
|
||||
ax_heat = fig.add_subplot(gs[0,0])
|
||||
ax_long = fig.add_subplot(gs[0,1], sharey=ax_heat)
|
||||
ax_trans = fig.add_subplot(gs[1,0], sharex=ax_heat)
|
||||
fig = plt.figure()
|
||||
gs = GridSpec(2, 2, width_ratios=[4, 1], height_ratios=[4, 1],
|
||||
hspace=0.3, wspace=0.1)
|
||||
|
||||
ax_heat = fig.add_subplot(gs[0, 0])
|
||||
ax_long = fig.add_subplot(gs[0, 1], sharey=ax_heat)
|
||||
ax_trans = fig.add_subplot(gs[1, 0], sharex=ax_heat)
|
||||
|
||||
# --------------------------
|
||||
# Heatmap Konturen
|
||||
# 2D-Konturen
|
||||
# --------------------------
|
||||
for E, color in zip(energies_MeV, colors):
|
||||
t = Z / X0_cm
|
||||
long_prof = longitudinal_profile(t, E, Ec_MeV, b) / X0_cm
|
||||
sigma_r = Rm_cm * (1 + 0.03 * t)
|
||||
trans_prof = np.exp(-(R**2)/(2*sigma_r**2))
|
||||
shower_2D = long_prof * trans_prof
|
||||
shower_norm = shower_2D / np.max(shower_2D) * 100
|
||||
for pct, ls in linestyles.items():
|
||||
ax_heat.contour(R, Z, shower_norm, levels=[pct], colors=[color], linestyles=[ls], linewidths=2)
|
||||
long_prof = dEdz(z, E)
|
||||
shower_2D = np.zeros((len(z), len(r)))
|
||||
rho_rz = rho_r(r)
|
||||
rho_rz[0] = rho_rz[1]
|
||||
|
||||
ax_heat.grid(True, linestyle='--', linewidth=0.5)
|
||||
ax_heat.set_ylabel("Depth z in cm", fontsize=fs)
|
||||
for i in range(len(z)):
|
||||
shower_2D[i, :] = long_prof[i] * rho_rz
|
||||
|
||||
levels = np.array(heatmap_levels)/100 * np.max(shower_2D)
|
||||
ax_heat.contour(R, Z, shower_2D,
|
||||
levels=levels,
|
||||
colors=[color]*len(levels),
|
||||
linestyles=linestyles)
|
||||
|
||||
# Achsen
|
||||
ax_heat.set_ylabel("Depth z in cm")
|
||||
ax_heat.set_ylim(max_depth_cm, 0)
|
||||
ax_heat.set_xlim(0,8)
|
||||
ax_heat.set_title("2D contourlines", fontsize=fs)
|
||||
ax_heat.tick_params(labelsize=fs)
|
||||
ax_heat.set_xlim(0, r_max)
|
||||
ax_heat.set_title("2D Energy density contours (normalized)")
|
||||
|
||||
# Molière Linien (Heatmap)
|
||||
molier_radii = [Rm_cm, 2*Rm_cm]
|
||||
molier_widths = [2, 2]
|
||||
for r_val, lw in zip(molier_radii, molier_widths):
|
||||
ax_heat.axvline(r_val, color='k', linestyle='-', linewidth=lw)
|
||||
ax_heat.axvline(-r_val, color='k', linestyle='-', linewidth=lw)
|
||||
# Molière-Linie
|
||||
molier_radii = [Rm_cm]
|
||||
for r_val in molier_radii:
|
||||
ax_heat.axvline(r_val, color='k', linestyle='-')
|
||||
|
||||
# --------------------------
|
||||
# Longitudinale Profile
|
||||
# Longitudinalprofile
|
||||
# --------------------------
|
||||
t = np.linspace(0, max_depth_cm/X0_cm, 400)
|
||||
t_plot = np.linspace(0, max_depth_cm / X0_cm, 400)
|
||||
for E, color in zip(energies_MeV, colors):
|
||||
profile = longitudinal_profile(t, E, Ec_MeV, b)/X0_cm
|
||||
ax_long.plot(profile, t*X0_cm, color=color, linewidth=2, label=f"{E/1000:.1f} GeV")
|
||||
profile = longitudinal_profile(t_plot, E, Ec_MeV, b) / X0_cm
|
||||
ax_long.plot(profile, t_plot*X0_cm, color=color, label=f"{E/1000:.1f} GeV")
|
||||
|
||||
ax_long.set_title("Longitudinal Profiles", fontsize=fs)
|
||||
ax_long.set_xlabel(r"$\mathrm{d}E/\mathrm{d}z$ in MeV/cm", fontsize=fs)
|
||||
ax_long.grid(True)
|
||||
ax_long.set_title("Longitudinal Profiles")
|
||||
ax_long.set_xlabel(r"$\mathrm{d}E/\mathrm{d}z$ in MeV/cm")
|
||||
ax_long.set_ylim(max_depth_cm, 0)
|
||||
ax_long.tick_params(labelsize=fs)
|
||||
ax_long.set_xlim(left=0)
|
||||
|
||||
# --------------------------
|
||||
# Transversale Profile (integrated over depth)
|
||||
# Transversales Profil
|
||||
# --------------------------
|
||||
n_bins = 80
|
||||
r_edges = np.linspace(0, r_max, n_bins+1)
|
||||
r_centers = 0.5*(r_edges[:-1]+r_edges[1:])
|
||||
z_grid = np.linspace(0, max_depth_cm, 400)
|
||||
dz = z_grid[1]-z_grid[0]
|
||||
|
||||
energy_lines = []
|
||||
dz = z[1] - z[0]
|
||||
for E, col in zip(energies_MeV, colors):
|
||||
prof_z = dEdz(z, E) # MeV/cm
|
||||
radial_profile = np.sum(prof_z[:, None] * rho_r(r)[None, :] * dz, axis=0) # MeV/cm²
|
||||
ln, = ax_trans.step(r, radial_profile, where='mid', color=col, linewidth=2)
|
||||
prof_z = dEdz(z_grid, E)
|
||||
rho = rho_r(r_centers)
|
||||
profile = np.sum(prof_z[:, None] * rho[None, :] * dz, axis=0)
|
||||
ln, = ax_trans.plot(r_centers, profile, color=col, label=f"{E/1000:.1f} GeV")
|
||||
energy_lines.append(ln)
|
||||
|
||||
# Molière Linien
|
||||
molier_lines = []
|
||||
for r_val, lw in zip(molier_radii, [2,2]):
|
||||
ln = ax_trans.axvline(r_val, color='k', linestyle='-', linewidth=lw)
|
||||
molier_lines.append(ln)
|
||||
ax_trans.axvline(-r_val, color='k', linestyle='-', linewidth=lw)
|
||||
for r_val in molier_radii:
|
||||
ax_trans.axvline(r_val, color='k', linestyle='-')
|
||||
|
||||
ax_trans.set_xlabel("Radius r in cm", fontsize=fs)
|
||||
ax_trans.set_ylabel(
|
||||
r"$\int \langle \mathrm{d}E/(\mathrm{d}z\,\mathrm{d}A) \rangle\,\mathrm{d}z$" + "\n in MeV/cm$^2$",
|
||||
fontsize=fs
|
||||
)
|
||||
ax_trans.set_title("Transverse energy deposition integrated over 10 cm BGO", fontsize=fs)
|
||||
ax_trans.grid(True)
|
||||
ax_trans.set_xlim(0,8)
|
||||
ax_trans.tick_params(labelsize=fs)
|
||||
|
||||
# --- Legenden ---
|
||||
leg1 = ax_heat.legend(energy_lines, [f"{E/1000:.1f} GeV" for E in energies_MeV],
|
||||
title="Initial Energy", fontsize=fs-2, title_fontsize=fs,
|
||||
loc="upper right", frameon=True, framealpha=0.85, facecolor="white")
|
||||
|
||||
leg2 = ax_trans.legend(molier_lines, [f"{r_val:.2f} cm" for r_val in molier_radii],
|
||||
title="Molière Radii", fontsize=fs-2, title_fontsize=fs,
|
||||
loc="upper right", frameon=True, framealpha=0.85, facecolor="white")
|
||||
|
||||
ax_heat.add_artist(leg1)
|
||||
ax_trans.add_artist(leg2)
|
||||
ax_trans.set_yscale("log")
|
||||
ax_trans.set_ylim(0.5, None)
|
||||
ax_trans.set_title("Transverse profiles integrated over 10cm")
|
||||
ax_trans.set_xlabel("Radius r in cm")
|
||||
ax_trans.set_ylabel(r"$\int \langle \mathrm{d}E/(\mathrm{d}z\,\mathrm{d}A) \rangle\, \mathrm{d}z$")
|
||||
|
||||
# --------------------------
|
||||
# Heatmap Legende (Contour)
|
||||
# Legenden
|
||||
# --------------------------
|
||||
line_legend = [Line2D([0],[0], color='k', linestyle=ls, lw=2) for ls in linestyles.values()]
|
||||
ax_heat.legend(line_legend, [f"{pct}\%" for pct in linestyles.keys()], title="Contour % of max",
|
||||
fontsize=fs-2, title_fontsize=fs, loc='lower right', frameon=True)
|
||||
|
||||
plt.tight_layout(rect=[0,0,0.95,0.95])
|
||||
plt.savefig("plots/shower_map_theory_depth.pdf")
|
||||
plt.savefig("plots/shower_map_theory_depth.png", dpi=300)
|
||||
ax_heat.legend(energy_lines, [f"{E/1000:.1f} GeV" for E in energies_MeV], title="Initial Energy", loc='upper right')
|
||||
plt.tight_layout()
|
||||
plotstyle.savefig("shower_map_theory_depth", category="BGO")
|
||||
plt.show()
|
||||
|
|
|
|||
|
|
@ -39,10 +39,18 @@ def dEdz(z, E):
|
|||
return dEdt / X0 # MeV/cm/nuc
|
||||
|
||||
# --------------------------
|
||||
# Transverse energy density
|
||||
# NKG-like transverse profile; safe at r=0
|
||||
# --------------------------
|
||||
s = 1.0
|
||||
def rho_r(r):
|
||||
return np.exp(-r**2/(2*Rm**2)) / (2*np.pi*Rm**2)
|
||||
x = r / Rm
|
||||
core = np.zeros_like(r)
|
||||
mask = r > 0
|
||||
core[mask] = (x[mask]**(s-2)) * ((1 + x[mask])**(s-4.5))
|
||||
if r[0] == 0:
|
||||
core[0] = ((1e-9 / Rm)**(s-2)) * ((1 + 1e-9 / Rm)**(s-4.5))
|
||||
norm = np.trapz(2*np.pi*r*core, r)
|
||||
return core / norm
|
||||
|
||||
# ============================================================
|
||||
# Plot 1: Longitudinal profile
|
||||
|
|
@ -64,15 +72,8 @@ plt.title("Longitudinal energy deposition in BGO (theory)", fontsize=fs+2)
|
|||
plt.xlim(0, z.max())
|
||||
plt.ylim(0, None)
|
||||
|
||||
leg = plt.legend(
|
||||
handles=lines,
|
||||
title="Initial Energy",
|
||||
fontsize=fs-1,
|
||||
title_fontsize=fs,
|
||||
frameon=True,
|
||||
framealpha=0.85,
|
||||
facecolor="white"
|
||||
)
|
||||
plt.legend(handles=lines, title="Initial Energy", fontsize=fs-1, title_fontsize=fs,
|
||||
frameon=True, framealpha=0.85, facecolor="white")
|
||||
|
||||
plt.grid(True)
|
||||
plt.xticks(fontsize=fs)
|
||||
|
|
@ -83,7 +84,7 @@ plt.savefig("plots/BGO_longitudinal_profile_normed.pdf")
|
|||
# ============================================================
|
||||
# Plot 2: Transverse profile at shower maximum
|
||||
# ============================================================
|
||||
r = np.linspace(0, 8, 300)
|
||||
r = np.linspace(0, 2.5, 300) # nur bis 2.5 cm
|
||||
|
||||
plt.figure(figsize=(12, 6))
|
||||
plt.style.use(['science'])
|
||||
|
|
@ -102,14 +103,14 @@ for E, lab, col in zip(energies, labels, colors):
|
|||
|
||||
# Geometry markers
|
||||
rm_line = plt.axvline(Rm, color='k', linestyle='--', linewidth=1.6, label=r"$R_M$")
|
||||
rm2_line = plt.axvline(2*Rm, color='k', linestyle=':', linewidth=1.6, label=r"$2R_M$")
|
||||
|
||||
plt.xlabel("Radius r in cm", fontsize=fs)
|
||||
plt.ylabel(r"$\langle dE/(\mathrm{d}z\,\mathrm{d}A) \rangle$" + "\nin MeV/cm$^2$", fontsize=fs)
|
||||
plt.title("Transverse energy density at shower maximum (theory)", fontsize=fs+2)
|
||||
|
||||
plt.xlim(0, r.max())
|
||||
plt.ylim(0, None)
|
||||
plt.xlim(0, 2.5)
|
||||
plt.yscale("log") # logarithmische y-Achse
|
||||
plt.ylim(0.1, 100000) # minimale y für log
|
||||
|
||||
# --- Legend 1: Energies ---
|
||||
leg1 = plt.legend(
|
||||
|
|
@ -125,21 +126,20 @@ leg1 = plt.legend(
|
|||
|
||||
# --- Legend 2: Geometry ---
|
||||
leg2 = plt.legend(
|
||||
handles=[rm_line, rm2_line],
|
||||
handles=[rm_line],
|
||||
fontsize=fs-1,
|
||||
loc="lower right",
|
||||
frameon=True,
|
||||
framealpha=0.85,
|
||||
facecolor="white"
|
||||
)
|
||||
|
||||
plt.gca().add_artist(leg1)
|
||||
|
||||
plt.grid(True)
|
||||
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
plt.tight_layout()
|
||||
plt.savefig("plots/BGO_transverse_profile_normed.pdf")
|
||||
plt.savefig("plots/BGO_transverse_profile_normed_log.pdf")
|
||||
|
||||
# ============================================================
|
||||
# Plot 3: Transverse profile integrated over entire BGO
|
||||
|
|
@ -159,15 +159,14 @@ for E, lab, col in zip(energies, labels, colors):
|
|||
|
||||
# Geometry markers
|
||||
rm_line = plt.axvline(Rm, color='k', linestyle='--', linewidth=1.6, label=r"$R_M$")
|
||||
rm2_line = plt.axvline(2*Rm, color='k', linestyle=':', linewidth=1.6, label=r"$2R_M$")
|
||||
|
||||
plt.xlabel("Radius r in cm", fontsize=fs)
|
||||
plt.ylabel(r"$\int \langle dE/(\mathrm{d}z\,\mathrm{d}A) \rangle\,dz$" + "\nin MeV/cm$^2$",fontsize=fs)
|
||||
plt.title("Transverse energy deposition integrated over 10 cm BGO (theory)", fontsize=fs+2)
|
||||
|
||||
|
||||
plt.xlim(0, r.max())
|
||||
plt.ylim(0, None)
|
||||
plt.xlim(0, 2.5)
|
||||
plt.yscale("log")
|
||||
plt.ylim(0.1, 100000)
|
||||
|
||||
# --- Legend 1: Energies ---
|
||||
leg1 = plt.legend(
|
||||
|
|
@ -183,20 +182,19 @@ leg1 = plt.legend(
|
|||
|
||||
# --- Legend 2: Geometry ---
|
||||
leg2 = plt.legend(
|
||||
handles=[rm_line, rm2_line],
|
||||
handles=[rm_line],
|
||||
fontsize=fs-1,
|
||||
loc="lower right",
|
||||
frameon=True,
|
||||
framealpha=0.85,
|
||||
facecolor="white"
|
||||
)
|
||||
|
||||
plt.gca().add_artist(leg1)
|
||||
|
||||
plt.grid(True)
|
||||
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
plt.tight_layout()
|
||||
plt.savefig("plots/BGO_transverse_profile_depth.pdf")
|
||||
plt.savefig("plots/BGO_transverse_profile_depth_log.pdf")
|
||||
|
||||
plt.show()
|
||||
|
|
|
|||
|
|
@ -28,7 +28,7 @@ colors = ["C0", "C1", "C2", "C3", "C4"]
|
|||
# -------------------------------------------------
|
||||
# Parameters
|
||||
# -------------------------------------------------
|
||||
r_max_cm = 8.0
|
||||
r_max_cm = 2.5
|
||||
n_bins = 100
|
||||
fs = 26
|
||||
|
||||
|
|
@ -84,54 +84,3 @@ plt.legend(title="Initial Energy", fontsize=fs-4, title_fontsize=fs-2, frameon=T
|
|||
plt.tight_layout()
|
||||
plt.savefig("plots/G4_transverse_profile_integrated.png", dpi=300)
|
||||
plt.show()
|
||||
|
||||
|
||||
# # --------- Dateien ---------
|
||||
# files = [
|
||||
# "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_100MeV_0.hits",
|
||||
# "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_200MeV_0.hits",
|
||||
# "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_500MeV_0.hits",
|
||||
# "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_1GeV_0.hits",
|
||||
# "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_2GeV_0.hits",
|
||||
# ]
|
||||
|
||||
# labels = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
|
||||
# colors = ["C0", "C1", "C2", "C3", "C4"]
|
||||
|
||||
# fs=18
|
||||
|
||||
# # --------- Parameter ---------
|
||||
# r_max_cm = 8.0
|
||||
# n_bins = 100
|
||||
|
||||
# # --------- Plot ---------
|
||||
# plt.figure(figsize=(12,6))
|
||||
|
||||
# for file, label, color in zip(files, labels, colors):
|
||||
# df = pd.read_csv(file, sep="\t")
|
||||
# df = df[df["r"] <= r_max_cm]
|
||||
|
||||
# r_edges = np.linspace(0, r_max_cm, n_bins + 1)
|
||||
# dr = r_edges[1] - r_edges[0]
|
||||
|
||||
# E_sum, _ = np.histogram(df["r"], bins=r_edges, weights=df["edep"])
|
||||
# counts, _ = np.histogram(df["r"], bins=r_edges)
|
||||
|
||||
# profile = np.zeros_like(E_sum)
|
||||
# mask = counts > 0
|
||||
# profile[mask] = E_sum[mask] / counts[mask] / dr # dE/dr [MeV/cm]
|
||||
|
||||
# r_centers = (r_edges[:-1] + r_edges[1:]) / 2
|
||||
# plt.step(r_centers, profile, lw=2, color=color, label=label)
|
||||
|
||||
# plt.xlabel("Radius r [cm]",fontsize=fs)
|
||||
# plt.ylabel("dE/dr [MeV/cm]",fontsize=fs)
|
||||
# plt.title("Transversal shower profile", fontsize=fs+2)
|
||||
# plt.xlim(0, r_max_cm)
|
||||
# plt.tick_params(axis='both', which='major', labelsize=fs-5)
|
||||
# plt.grid(True, ls="--", lw=0.5)
|
||||
# plt.legend(title="Initial Energy")
|
||||
|
||||
# plt.tight_layout()
|
||||
# plt.savefig("plots/shower_transverse.png", dpi=300)
|
||||
# plt.show()
|
||||
|
|
|
|||
180
shower_trans_comp_slices.py
Normal file
180
shower_trans_comp_slices.py
Normal file
|
|
@ -0,0 +1,180 @@
|
|||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Transverse shower profiles in BGO
|
||||
Comparison: Theory (NKG, dashed) vs Geant4 simulation (step)
|
||||
|
||||
Plot:
|
||||
- 5 slices (0-2,2-4,4-6,6-8,8-10 cm)
|
||||
- 1 sum plot (0-10 cm)
|
||||
- Legend: Initial energy + Theory vs Geant4
|
||||
- Slice-wise Integration für Theorie
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
from scipy.stats import gamma
|
||||
import scienceplots
|
||||
from matplotlib.lines import Line2D
|
||||
|
||||
# ============================================================
|
||||
# Style & Font
|
||||
# ============================================================
|
||||
plt.style.use(['science'])
|
||||
fs = 26
|
||||
|
||||
# ============================================================
|
||||
# Material (BGO)
|
||||
# ============================================================
|
||||
X0 = 1.118
|
||||
Ec = 10.5
|
||||
b = 0.5
|
||||
Rm = 2.26
|
||||
|
||||
# ============================================================
|
||||
# Inputs
|
||||
# ============================================================
|
||||
files = [
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_100MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_200MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_500MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_1GeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_2GeV_0.hits",
|
||||
]
|
||||
|
||||
energies = [100, 200, 500, 1000, 2000]
|
||||
labels = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
|
||||
colors = ["C0", "C1", "C2", "C3", "C4"]
|
||||
|
||||
# ============================================================
|
||||
# Theory functions
|
||||
# ============================================================
|
||||
def tmax(E):
|
||||
return np.log(E / Ec) - 0.5
|
||||
|
||||
def dEdz(z, E):
|
||||
t = z / X0
|
||||
a = 1 + b * tmax(E)
|
||||
return E * gamma.pdf(t, a, scale=1/b) # MeV/cm
|
||||
|
||||
# NKG transverse profile
|
||||
s = 1.0
|
||||
def rho_r(r):
|
||||
x = r / Rm
|
||||
core = np.zeros_like(r)
|
||||
mask = r > 0
|
||||
core[mask] = (x[mask]**(s-2)) * ((1 + x[mask])**(s-4.5))
|
||||
if r[0] == 0:
|
||||
core[0] = ((1e-9 / Rm)**(s-2)) * ((1 + 1e-9 / Rm)**(s-4.5))
|
||||
norm = np.trapz(2*np.pi*r*core, r)
|
||||
return core / norm
|
||||
|
||||
# ============================================================
|
||||
# Radial binning
|
||||
# ============================================================
|
||||
r_max_cm = 2.5
|
||||
n_bins = 80
|
||||
r_edges = np.linspace(0, r_max_cm, n_bins+1)
|
||||
r_centers = 0.5*(r_edges[:-1]+r_edges[1:])
|
||||
A_ring = np.pi * (r_edges[1:]**2 - r_edges[:-1]**2)
|
||||
|
||||
# ============================================================
|
||||
# Slices
|
||||
# ============================================================
|
||||
z_slices = [(0,2),(2,4),(4,6),(6,8),(8,10)]
|
||||
slice_labels = ["0-2 cm","2-4 cm","4-6 cm","6-8 cm","8-10 cm"]
|
||||
|
||||
# ============================================================
|
||||
# Figure
|
||||
# ============================================================
|
||||
fig, axes = plt.subplots(6,1, figsize=(13,16), sharex=True,
|
||||
gridspec_kw={"height_ratios":[1,1,1,1,1,1],"hspace":0.05})
|
||||
|
||||
z_grid = np.linspace(0,10,400)
|
||||
dz = z_grid[1]-z_grid[0]
|
||||
rho = rho_r(r_centers)
|
||||
|
||||
# Handle-Listen für Energy-Legende (oben)
|
||||
energy_lines = []
|
||||
|
||||
for file, E, lab, col in zip(files, energies, labels, colors):
|
||||
df = pd.read_csv(file, sep="\t")
|
||||
df = df[df["r"] <= r_max_cm]
|
||||
N_events = df["event"].nunique()
|
||||
|
||||
# Für Startenergie-Legende
|
||||
energy_lines.append(Line2D([0],[0], color=col, lw=2))
|
||||
|
||||
# --- Slices ---
|
||||
for i, (z1, z2) in enumerate(z_slices):
|
||||
ax = axes[i]
|
||||
|
||||
# Simulation
|
||||
df_slice = df[(df["z"]>=z1)&(df["z"]<z2)]
|
||||
E_sum_slice,_ = np.histogram(df_slice["r"], bins=r_edges, weights=df_slice["edep"])
|
||||
sim_profile = E_sum_slice / (N_events*A_ring)
|
||||
|
||||
# Theorie Slice-wise
|
||||
mask = (z_grid >= z1) & (z_grid < z2)
|
||||
theory_profile = np.sum(dEdz(z_grid[mask], E)[:, None] * rho[None, :] * dz, axis=0)
|
||||
|
||||
ax.plot(r_centers, theory_profile, '--', lw=2, color=col)
|
||||
ax.step(r_centers, sim_profile, where='mid', lw=2, color=col)
|
||||
|
||||
ax.set_yscale('log')
|
||||
ax.set_ylim(0.05, 30000)
|
||||
ax.set_ylabel(slice_labels[i], fontsize=fs-6)
|
||||
ax.grid(True, which='both', ls='--', lw=0.6)
|
||||
|
||||
# --- Summe 0-10 cm ---
|
||||
ax = axes[5]
|
||||
# Simulation
|
||||
E_sum, _ = np.histogram(df["r"], bins=r_edges, weights=df["edep"])
|
||||
sim_profile = E_sum / (N_events * A_ring)
|
||||
# Theorie über Gesamt 0-10 cm
|
||||
theory_profile = np.sum(dEdz(z_grid, E)[:, None] * rho[None, :] * dz, axis=0)
|
||||
|
||||
ax.plot(r_centers, theory_profile, '--', lw=2, color=col)
|
||||
ax.step(r_centers, sim_profile, where='mid', lw=2, color=col)
|
||||
|
||||
ax.set_yscale('log')
|
||||
ax.set_ylim(0.05, 30000)
|
||||
ax.set_ylabel('Sum', fontsize=fs-6)
|
||||
ax.grid(True, which='both', ls='--', lw=0.6)
|
||||
|
||||
# ============================================================
|
||||
# Cosmetics
|
||||
# ============================================================
|
||||
axes[-1].set_xlabel('Radius r in cm', fontsize=fs)
|
||||
for ax in axes:
|
||||
ax.set_xlim(0,r_max_cm)
|
||||
ax.tick_params(labelsize=fs-2)
|
||||
|
||||
fig.text(0.04, 0.5,
|
||||
r"$\left\langle \int \frac{\mathrm{d}E}{\mathrm{d}z\,\mathrm{d}A} \, \mathrm{d}z \right\rangle$ in MeV/cm$^2$",
|
||||
va='center', ha='center', rotation='vertical', fontsize=fs)
|
||||
|
||||
# ============================================================
|
||||
# Überschrift
|
||||
# ============================================================
|
||||
fig.suptitle("Transverse shower profiles in BGO", fontsize=fs+2, y=0.94)
|
||||
|
||||
# ============================================================
|
||||
# Legend: Initialenergie zwischen Titel und oberstem Plot
|
||||
# ============================================================
|
||||
fig.legend(energy_lines, labels, loc='upper center', bbox_to_anchor=(0.5, 0.92),
|
||||
ncol=5, fontsize=fs-10, frameon=True)
|
||||
|
||||
# ============================================================
|
||||
# Legend: Theory vs Geant4 Simulation
|
||||
# ============================================================
|
||||
sim_lines = [
|
||||
Line2D([0],[0], color='k', lw=2, linestyle='-', label='Geant4 simulation'),
|
||||
Line2D([0],[0], color='k', lw=2, linestyle='--', label='Theory')
|
||||
]
|
||||
axes[0].legend(handles=sim_lines, loc='upper right', ncol=2, fontsize=fs-10, frameon=True)
|
||||
|
||||
plt.tight_layout(rect=[0,0,1,0.85])
|
||||
plt.savefig("plots/shower_trans_comp_slices.pdf", dpi=300)
|
||||
plt.show()
|
||||
224
shower_trans_comparison.py
Normal file
224
shower_trans_comparison.py
Normal file
|
|
@ -0,0 +1,224 @@
|
|||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Transverse shower profiles in BGO
|
||||
Comparison: Theory (NKG, dashed) vs Geant4 simulation (step)
|
||||
Residuals shown logarithmically
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
from scipy.stats import gamma
|
||||
from matplotlib.lines import Line2D
|
||||
import sys
|
||||
from pathlib import Path
|
||||
|
||||
sys.path.append(str(Path(__file__).resolve().parents[1]))
|
||||
import plotstyle
|
||||
|
||||
import scienceplots
|
||||
|
||||
# ============================================================
|
||||
# Style & Font
|
||||
# ============================================================
|
||||
plotstyle.set_style("paper") # nutzt Slides-Mode
|
||||
fs = plt.rcParams["font.size"]
|
||||
|
||||
# ============================================================
|
||||
# Material (BGO)
|
||||
# ============================================================
|
||||
X0 = 1.118
|
||||
Ec = 10.5
|
||||
b = 0.5
|
||||
Rm = 2.26
|
||||
|
||||
# ============================================================
|
||||
# Inputs
|
||||
# ============================================================
|
||||
files = [
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_100MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_200MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_500MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_1GeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_2GeV_0.hits",
|
||||
]
|
||||
|
||||
energies = [100, 200, 500, 1000, 2000]
|
||||
labels = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
|
||||
colors = ["C0", "C1", "C2", "C3", "C4"]
|
||||
|
||||
# ============================================================
|
||||
# Theory functions
|
||||
# ============================================================
|
||||
def tmax(E):
|
||||
return np.log(E / Ec) - 0.5
|
||||
|
||||
def dEdz(z, E):
|
||||
t = z / X0
|
||||
a = 1 + b * tmax(E)
|
||||
return E * gamma.pdf(t, a, scale=1/b) / X0 # MeV/cm
|
||||
|
||||
# NKG transverse profile
|
||||
s = 1.0
|
||||
def rho_r(r):
|
||||
x = r / Rm
|
||||
core = np.zeros_like(r)
|
||||
mask = r > 0
|
||||
core[mask] = (x[mask]**(s-2)) * ((1 + x[mask])**(s-4.5))
|
||||
if r[0] == 0:
|
||||
core[0] = ((1e-9 / Rm)**(s-2)) * ((1 + 1e-9 / Rm)**(s-4.5))
|
||||
norm = np.trapz(2*np.pi*r*core, r)
|
||||
return core / norm
|
||||
|
||||
# ============================================================
|
||||
# Common binning
|
||||
# ============================================================
|
||||
r_max_cm = 2.5
|
||||
n_bins = 80
|
||||
|
||||
r_edges = np.linspace(0, r_max_cm, n_bins + 1)
|
||||
r_centers = 0.5 * (r_edges[:-1] + r_edges[1:])
|
||||
r_in = r_edges[:-1]
|
||||
r_out = r_edges[1:]
|
||||
A_ring = np.pi * (r_out**2 - r_in**2)
|
||||
|
||||
# ============================================================
|
||||
# Legend styles
|
||||
# ============================================================
|
||||
style_legend = [
|
||||
Line2D([0], [0], color="k", lw=2, linestyle="--", label="Theory"),
|
||||
Line2D([0], [0], color="k", lw=2, linestyle="-", label="Geant4 simulation"),
|
||||
]
|
||||
|
||||
# ============================================================
|
||||
# FIGURE 1 — Shower Maximum
|
||||
# ============================================================
|
||||
fig, (ax, ax_res) = plt.subplots(
|
||||
2, 1,
|
||||
sharex=True,
|
||||
gridspec_kw={"height_ratios": [3, 1], "hspace": 0.08}
|
||||
)
|
||||
|
||||
energy_lines = []
|
||||
dz_slice = 0.5
|
||||
|
||||
for file, E, lab, col in zip(files, energies, labels, colors):
|
||||
|
||||
df = pd.read_csv(file, sep="\t")
|
||||
N_events = df["event"].nunique()
|
||||
|
||||
z0 = tmax(E) * X0
|
||||
|
||||
df_slice = df[
|
||||
(df["z"] >= z0 - dz_slice/2) &
|
||||
(df["z"] <= z0 + dz_slice/2) &
|
||||
(df["r"] <= r_max_cm)
|
||||
]
|
||||
|
||||
E_sum, _ = np.histogram(df_slice["r"], bins=r_edges, weights=df_slice["edep"])
|
||||
sim_profile = E_sum / (N_events * dz_slice * A_ring) # MeV/cm^2
|
||||
|
||||
# Theory
|
||||
theory_profile = dEdz(z0, E) * rho_r(r_centers) # MeV/cm^2
|
||||
|
||||
ax.plot(r_centers, theory_profile, linestyle="--", linewidth=2, color=col)
|
||||
ax.step(r_centers, sim_profile, where="mid", linewidth=2, color=col)
|
||||
|
||||
energy_lines.append(Line2D([0], [0], color=col, lw=2, linestyle="-"))
|
||||
|
||||
# Residuals
|
||||
ax_res.step(r_centers, (sim_profile - theory_profile) / theory_profile, color=col, linewidth=1.5)
|
||||
|
||||
ax.set_yscale("log")
|
||||
ax.set_xlim(0, r_max_cm)
|
||||
ax.set_ylabel(
|
||||
r"$\left\langle \frac{\mathrm{d}E}{\mathrm{d}z\,\mathrm{d}A} \right\rangle$ in MeV/cm$^2$",
|
||||
fontsize=fs
|
||||
)
|
||||
ax.set_title("Transverse shower profile at shower maximum (BGO)")
|
||||
ax_res.set_xlim(0, r_max_cm)
|
||||
ax_res.set_xlabel("Radius r in cm")
|
||||
ax_res.set_ylabel(r"$\frac{\mathrm{G4sim} - \mathrm{Theory}}{\mathrm{Theory}}$",)
|
||||
ax_res.axhline(0, color='k', linestyle=':', linewidth=1.5)
|
||||
|
||||
# Legends
|
||||
leg1 = ax.legend(
|
||||
energy_lines, labels,
|
||||
title="Initial Energy",
|
||||
loc="upper right"
|
||||
)
|
||||
|
||||
leg2 = ax.legend(
|
||||
handles=style_legend,
|
||||
loc="lower left",
|
||||
)
|
||||
ax.add_artist(leg1)
|
||||
|
||||
plotstyle.savefig("shower_trans_comparison_max", category="BGO")
|
||||
|
||||
# ============================================================
|
||||
# FIGURE 2 — Integrated over crystal
|
||||
# ============================================================
|
||||
fig, (ax, ax_res) = plt.subplots(
|
||||
2, 1,
|
||||
sharex=True,
|
||||
gridspec_kw={"height_ratios": [3, 1], "hspace": 0.05}
|
||||
)
|
||||
|
||||
z_grid = np.linspace(0, 10, 400)
|
||||
dz = z_grid[1] - z_grid[0]
|
||||
|
||||
for file, E, lab, col in zip(files, energies, labels, colors):
|
||||
|
||||
df = pd.read_csv(file, sep="\t")
|
||||
df = df[df["r"] <= r_max_cm]
|
||||
N_events = df["event"].nunique()
|
||||
|
||||
# Simulation
|
||||
E_sum, _ = np.histogram(df["r"], bins=r_edges, weights=df["edep"])
|
||||
sim_profile = E_sum / (N_events * A_ring) # MeV/cm^2
|
||||
|
||||
# Theory
|
||||
prof_z = dEdz(z_grid, E)
|
||||
theory_profile = np.sum(
|
||||
prof_z[:, None] * rho_r(r_centers)[None, :] * dz,
|
||||
axis=0
|
||||
)
|
||||
|
||||
ax.plot(r_centers, theory_profile, linestyle="--", linewidth=2, color=col)
|
||||
ax.step(r_centers, sim_profile, where="mid", linewidth=2, color=col)
|
||||
ax_res.step(r_centers, (sim_profile - theory_profile)/theory_profile, color=col, linewidth=1.5)
|
||||
|
||||
ax.set_yscale("log")
|
||||
ax.set_ylim(0.05, 30000)
|
||||
ax.set_xlim(0, r_max_cm)
|
||||
ax.set_ylabel(
|
||||
r"$\left\langle \int \frac{\mathrm{d}E}{\mathrm{d}z\,\mathrm{d}A} \, \mathrm{d}z \right\rangle$ in MeV/cm$^2$",
|
||||
)
|
||||
ax.set_title("Transverse shower profile integrated over 10 cm BGO")
|
||||
|
||||
ax_res.set_xlim(0, r_max_cm)
|
||||
ax_res.set_xlabel("Radius r in cm")
|
||||
ax_res.set_ylabel(
|
||||
r"$\frac{\mathrm{G4sim} - \mathrm{Theory}}{\mathrm{Theory}}$",
|
||||
)
|
||||
ax_res.axhline(0, color='k', linestyle=':', linewidth=1.5)
|
||||
|
||||
# Legends
|
||||
leg1 = ax.legend(
|
||||
energy_lines, labels,
|
||||
title="Initial Energy",
|
||||
title_fontsize=fs-2,
|
||||
loc="upper right",
|
||||
)
|
||||
|
||||
leg2 = ax.legend(
|
||||
handles=style_legend,
|
||||
loc="lower left",
|
||||
)
|
||||
ax.add_artist(leg1)
|
||||
|
||||
plotstyle.savefig("shower_trans_comparison_integrated", category="BGO")
|
||||
|
||||
plt.show()
|
||||
|
|
@ -39,7 +39,7 @@ colors = ["C0", "C1", "C2", "C3", "C4"]
|
|||
# -------------------------------------------------
|
||||
# Analysis parameters
|
||||
# -------------------------------------------------
|
||||
r_max_cm = 8.0
|
||||
r_max_cm = 2.5
|
||||
n_bins = 80
|
||||
|
||||
dz_slice = 0.5 # cm (slice thickness around shower maximum)
|
||||
|
|
@ -111,11 +111,113 @@ plt.title("Transverse energy density at shower maximum (Geant4)", fontsize=fs+2)
|
|||
|
||||
plt.xlim(0, r_max_cm)
|
||||
plt.yscale("log")
|
||||
plt.grid(True, ls="--", lw=0.6)
|
||||
plt.grid(True)
|
||||
plt.tick_params(labelsize=fs-2)
|
||||
plt.legend(title="Initial Energy", fontsize=fs-4, title_fontsize=fs-2, frameon=True, framealpha=0.85, facecolor="white")
|
||||
|
||||
|
||||
plt.tight_layout()
|
||||
plt.savefig("plots/G4_transverse_profile_showermax.png", dpi=300)
|
||||
|
||||
# --------------------------
|
||||
# Fixed z positions
|
||||
# --------------------------
|
||||
z_positions = [2.0, 4.0, 6.0] # cm
|
||||
dz_slice = 0.5 # cm slice thickness
|
||||
n_slices = len(z_positions)
|
||||
|
||||
# --------------------------
|
||||
# Figure
|
||||
# --------------------------
|
||||
plt.style.use(['science'])
|
||||
fig, axes = plt.subplots(n_slices, 1, figsize=(10, 10), sharex=True)
|
||||
plt.subplots_adjust(hspace=0.12)
|
||||
|
||||
fig.suptitle(
|
||||
"Transverse energy density at fixed depths in BGO (Geant4)",
|
||||
fontsize=fs+2,
|
||||
y=0.95
|
||||
)
|
||||
|
||||
# Dummy lines for legend (energy only)
|
||||
dummy_lines = [axes[0].plot([], [], color=c, linewidth=2)[0] for c in colors]
|
||||
axes[0].legend(
|
||||
dummy_lines,
|
||||
[f"{E} MeV" for E in energies],
|
||||
ncol=5,
|
||||
fontsize=fs-10,
|
||||
frameon=True,
|
||||
framealpha=0.85,
|
||||
loc="upper center",
|
||||
bbox_to_anchor=(0.5, 1.25)
|
||||
)
|
||||
|
||||
# Shared y-label
|
||||
fig.text(
|
||||
0.02, 0.5,
|
||||
r"$\langle \mathrm{d}E/(\mathrm{d}z\,\mathrm{d}A) \rangle$ in MeV/cm$^3$",
|
||||
va="center",
|
||||
rotation="vertical",
|
||||
fontsize=fs
|
||||
)
|
||||
|
||||
# --------------------------
|
||||
# Loop over z slices
|
||||
# --------------------------
|
||||
for ax, z_ref in zip(axes, z_positions):
|
||||
|
||||
for file, E, color in zip(files, energies, colors):
|
||||
df = pd.read_csv(file, sep="\t")
|
||||
|
||||
# Number of primary events
|
||||
N_events = df["event"].nunique()
|
||||
|
||||
# z slice around fixed depth
|
||||
df_slice = df[
|
||||
(df["z"] >= z_ref - dz_slice/2) &
|
||||
(df["z"] <= z_ref + dz_slice/2)
|
||||
]
|
||||
|
||||
df_slice = df_slice[df_slice["r"] <= r_max_cm]
|
||||
|
||||
# Radial histogram
|
||||
r_edges = np.linspace(0.0, r_max_cm, n_bins + 1)
|
||||
r_centers = 0.5 * (r_edges[:-1] + r_edges[1:])
|
||||
|
||||
E_sum, _ = np.histogram(
|
||||
df_slice["r"],
|
||||
bins=r_edges,
|
||||
weights=df_slice["edep"]
|
||||
)
|
||||
|
||||
# Ring areas
|
||||
r_in = r_edges[:-1]
|
||||
r_out = r_edges[1:]
|
||||
A_ring = np.pi * (r_out**2 - r_in**2)
|
||||
|
||||
# Normalised profile
|
||||
profile = E_sum / (N_events * dz_slice * A_ring)
|
||||
|
||||
ax.step(
|
||||
r_centers,
|
||||
profile,
|
||||
where="mid",
|
||||
lw=2,
|
||||
color=color
|
||||
)
|
||||
|
||||
ax.set_xlim(0, r_max_cm)
|
||||
ax.set_yscale("log")
|
||||
ax.grid(True, ls="--", lw=0.6)
|
||||
ax.tick_params(labelsize=fs-2)
|
||||
|
||||
# z-label on the left of each panel
|
||||
ax.set_ylabel(f"z = {z_ref:.0f} cm", fontsize=fs)
|
||||
|
||||
# Bottom x-label
|
||||
axes[-1].set_xlabel("Radius r in cm", fontsize=fs)
|
||||
|
||||
plt.tight_layout(rect=[0.06, 0.03, 0.98, 0.93])
|
||||
plt.savefig("plots/G4_transverse_profile_fixed_z.png", dpi=300)
|
||||
plt.show()
|
||||
|
||||
|
|
|
|||
|
|
@ -9,6 +9,10 @@ import numpy as np
|
|||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
import scienceplots
|
||||
from matplotlib.lines import Line2D
|
||||
|
||||
plt.style.use(['science'])
|
||||
fs = 26
|
||||
|
||||
# --------------------------
|
||||
# Input files and parameters
|
||||
|
|
@ -30,7 +34,7 @@ fs = 24
|
|||
# Analysis parameters
|
||||
# --------------------------
|
||||
z_max_cm = 10.0
|
||||
r_max_cm = 8.0
|
||||
r_max_cm = 2.5
|
||||
dz_layer = 2.0 # cm
|
||||
n_bins_r = 100
|
||||
|
||||
|
|
@ -44,20 +48,8 @@ dr = r_edges[1] - r_edges[0]
|
|||
# --------------------------
|
||||
# Figure
|
||||
# --------------------------
|
||||
plt.style.use(['science'])
|
||||
fig, axes = plt.subplots(n_layers+1, 1, figsize=(10,12), sharex=True)
|
||||
plt.subplots_adjust(hspace=0.1)
|
||||
fig.suptitle("Transverse energy deposition in BGO layers (Geant4)", fontsize=fs+2, y=0.95)
|
||||
|
||||
# Dummy-Linien für Startenergie-Legende
|
||||
dummy_lines = [axes[0].plot([], [], color=c, linewidth=2)[0] for c in colors]
|
||||
axes[0].legend(dummy_lines, [f"{E} MeV" for E in energies],
|
||||
ncol=5, fontsize=fs-10, frameon=True, framealpha=0.85,
|
||||
loc='upper center', bbox_to_anchor=(0.5, 1.05))
|
||||
|
||||
# Gemeinsames Y-Label
|
||||
fig.text(0.02, 0.5, r"$\langle \mathrm{d}E/(\mathrm{d}z\,\mathrm{d}A) \rangle$ in MeV/cm$^2$",
|
||||
va='center', rotation='vertical', fontsize=fs)
|
||||
|
||||
# --------------------------
|
||||
# Layerwise plotting with integrals
|
||||
|
|
@ -96,8 +88,8 @@ for i, (z_min, z_max) in enumerate(layers + [(0,z_max_cm)]):
|
|||
|
||||
ax.set_xlim(0, r_max_cm)
|
||||
ax.set_yscale('log')
|
||||
#ax.set_ylim(0.0009, None)
|
||||
ax.grid(True)
|
||||
ax.set_ylim(0.01, 1e5)
|
||||
ax.grid(True, which='both', linestyle='--', linewidth=0.5)
|
||||
ax.tick_params(labelsize=fs-2)
|
||||
|
||||
if i < n_layers:
|
||||
|
|
@ -109,83 +101,26 @@ for i, (z_min, z_max) in enumerate(layers + [(0,z_max_cm)]):
|
|||
# Legende innerhalb der Achse
|
||||
ax.legend(fontsize=12, frameon=True, framealpha=0.85, facecolor="white")
|
||||
|
||||
# Cosmetics
|
||||
axes[-1].set_xlabel('Radius r in cm', fontsize=fs)
|
||||
for ax in axes:
|
||||
ax.set_xlim(0,r_max_cm)
|
||||
ax.tick_params(labelsize=fs-2)
|
||||
|
||||
# Gemeinsames y-Label
|
||||
fig.text(0.04, 0.5,
|
||||
r"$\left\langle \int \frac{dE}{dz\,dA} dz \right\rangle$ [MeV/cm$^2$]",
|
||||
va='center', ha='center', rotation='vertical', fontsize=fs)
|
||||
|
||||
# Überschrift
|
||||
fig.suptitle("Transverse shower profiles in BGO (Geant4 Simulation)", fontsize=fs+2, y=0.98)
|
||||
|
||||
# Startenergie-Legende oberhalb oberstem Plot
|
||||
energy_lines = [Line2D([0],[0], color=c,lw=2) for c in colors]
|
||||
fig.legend(energy_lines, labels, loc='upper center', bbox_to_anchor=(0.5,0.95),
|
||||
ncol=5, fontsize=fs-10, frameon=True)
|
||||
|
||||
plt.tight_layout(rect=[0.05,0.03,0.97,0.93])
|
||||
plt.savefig("plots/G4_transverse_layers_sum.pdf")
|
||||
plt.savefig("plots/G4_transverse_layers_sum.png", dpi=300)
|
||||
plt.show()
|
||||
|
||||
|
||||
|
||||
|
||||
# #!/usr/bin/env python3
|
||||
# # -*- coding: utf-8 -*-
|
||||
# """
|
||||
# transversal shower profile form Genatz4 simulations slicewise
|
||||
# """
|
||||
|
||||
# import numpy as np
|
||||
# import pandas as pd
|
||||
# import matplotlib.pyplot as plt
|
||||
|
||||
# # --------- Dateien ---------
|
||||
# files = [
|
||||
# "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_100MeV_0.hits",
|
||||
# "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_200MeV_0.hits",
|
||||
# "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_500MeV_0.hits",
|
||||
# "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_1GeV_0.hits",
|
||||
# "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_2GeV_0.hits",
|
||||
# ]
|
||||
|
||||
# labels = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
|
||||
# colors = ["C0", "C1", "C2", "C3", "C4"]
|
||||
|
||||
# fs = 18 # Schriftgröße
|
||||
# r_max_cm = 8.0
|
||||
# n_bins = 100
|
||||
# n_particles = 10000 # Anzahl der simulierten Teilchen
|
||||
|
||||
# # --------- Subplots ---------
|
||||
# fig, axes = plt.subplots(5, 1, figsize=(8, 12), sharex=True)
|
||||
# z_slices = [(0, 2), (2, 4), (4, 6), (6, 8)]
|
||||
# dz_slice = 2.0 # Breite der Slice in cm
|
||||
|
||||
# # --------- Plotten ---------
|
||||
# for file, label, color in zip(files, labels, colors):
|
||||
# df = pd.read_csv(file, sep="\t")
|
||||
# df = df[df["r"] <= r_max_cm]
|
||||
|
||||
# r_edges = np.linspace(0, r_max_cm, n_bins + 1)
|
||||
# dr = r_edges[1] - r_edges[0]
|
||||
# r_centers = (r_edges[:-1] + r_edges[1:]) / 2
|
||||
|
||||
# # --------- Obere 4 Slices ---------
|
||||
# for i, (z_min, z_max) in enumerate(z_slices):
|
||||
# df_slice = df[(df["z"] >= z_min) & (df["z"] < z_max)]
|
||||
# E_sum, _ = np.histogram(df_slice["r"], bins=r_edges, weights=df_slice["edep"])
|
||||
# # Mittelwert pro Teilchen pro cm
|
||||
# profile = E_sum / (n_particles * dz_slice)
|
||||
# axes[i].step(r_centers, profile, lw=2, color=color, label=label if i==0 else "")
|
||||
# axes[i].set_ylim(0, None)
|
||||
# axes[i].tick_params(axis='both', which='major', labelsize=fs-4)
|
||||
# # kleines y-Label rechts mit Slice
|
||||
# axes[i].text(r_max_cm*1.01, axes[i].get_ylim()[1]*0.9, f"{z_min}-{z_max} cm", rotation=0, va="top", fontsize=fs-6)
|
||||
|
||||
# # --------- Unterer Plot: gesamte Summe 0-8cm ---------
|
||||
# df_all = df[df["r"] <= r_max_cm]
|
||||
# E_sum_total, _ = np.histogram(df_all["r"], bins=r_edges, weights=df_all["edep"])
|
||||
# profile_total = E_sum_total / (n_particles * 8.0) # mittlerer Energieverlust pro cm
|
||||
# axes[4].step(r_centers, profile_total, lw=2, color='k')
|
||||
# axes[4].set_ylim(0, None)
|
||||
# axes[4].tick_params(axis='both', which='major', labelsize=fs-4)
|
||||
# axes[4].text(r_max_cm*1.01, axes[4].get_ylim()[1]*0.9, "0-8 cm", rotation=0, va="top", fontsize=fs-6)
|
||||
|
||||
# # --------- Gemeinsames y-Label ---------
|
||||
# fig.text(0.02, 0.5, "Energy loss [MeV/cm]", va='center', rotation='vertical', fontsize=fs)
|
||||
|
||||
# # --------- Achsen & Legende ---------
|
||||
# axes[-1].set_xlabel("Radius r [cm]", fontsize=fs)
|
||||
# axes[0].legend(title="Initial Energy", loc="upper right", fontsize=fs-4)
|
||||
|
||||
# plt.tight_layout(rect=[0.05,0.03,1,0.97])
|
||||
# plt.savefig("plots/shower_transverse_slices.png", dpi=300)
|
||||
# plt.show()
|
||||
plt.show()
|
||||
121
shower_trans_slices2.py
Normal file
121
shower_trans_slices2.py
Normal file
|
|
@ -0,0 +1,121 @@
|
|||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Transverse shower profiles in BGO
|
||||
Nur Geant4 Simulation
|
||||
Legenden: obere Legende = Startenergien, rechts pro Subplot = integrierte Energie pro Slice
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
import scienceplots
|
||||
from matplotlib.lines import Line2D
|
||||
|
||||
plt.style.use(['science'])
|
||||
fs = 26
|
||||
|
||||
# Radial binning
|
||||
r_max_cm = 2.5
|
||||
n_bins = 80
|
||||
r_edges = np.linspace(0, r_max_cm, n_bins+1)
|
||||
r_centers = 0.5*(r_edges[:-1] + r_edges[1:])
|
||||
r_in = r_edges[:-1]
|
||||
r_out = r_edges[1:]
|
||||
A_ring = np.pi*(r_out**2 - r_in**2)
|
||||
|
||||
# Slices
|
||||
z_slices = [(0,2), (2,4), (4,6), (6,8), (8,10)]
|
||||
slice_labels = ["0-2 cm", "2-4 cm", "4-6 cm", "6-8 cm", "8-10 cm"]
|
||||
|
||||
# Inputs
|
||||
files = [
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_100MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_200MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_500MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_1GeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_2GeV_0.hits",
|
||||
]
|
||||
labels = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
|
||||
colors = ["C0", "C1", "C2", "C3", "C4"]
|
||||
|
||||
# Figure
|
||||
fig, axes = plt.subplots(6, 1, figsize=(13,16), sharex=True,
|
||||
gridspec_kw={"height_ratios":[1,1,1,1,1,1], "hspace":0.05})
|
||||
|
||||
# Bereite Listen für Handles und Labels pro Subplot vor (nur integrierte Energie)
|
||||
handles_slices = [[] for _ in range(len(z_slices)+1)]
|
||||
labels_slices = [[] for _ in range(len(z_slices)+1)]
|
||||
|
||||
# Handle-Liste für Startenergie-Legende oben
|
||||
energy_lines = []
|
||||
|
||||
for file, E, col, label in zip(files, [100,200,500,1000,2000], colors, labels):
|
||||
df = pd.read_csv(file, sep="\t")
|
||||
df = df[df["r"] <= r_max_cm]
|
||||
N_events = df["event"].nunique()
|
||||
|
||||
# Für Startenergie-Legende
|
||||
energy_lines.append(Line2D([0], [0], color=col, lw=2))
|
||||
|
||||
for i, (z1, z2) in enumerate(z_slices):
|
||||
ax = axes[i]
|
||||
df_slice = df[(df["z"] >= z1) & (df["z"] < z2)]
|
||||
E_sum_slice, _ = np.histogram(df_slice["r"], bins=r_edges, weights=df_slice["edep"])
|
||||
sim_profile = E_sum_slice / (N_events * A_ring)
|
||||
|
||||
# Physikalisch korrekte Slice-Energie
|
||||
E_slice_total = np.sum(E_sum_slice) / N_events
|
||||
E_slice_total = min(E_slice_total, E)
|
||||
|
||||
line = ax.step(r_centers, sim_profile, where='mid', lw=2, color=col)[0]
|
||||
|
||||
# Nur integrierte Energie als Subplot-Legende
|
||||
handles_slices[i].append(line)
|
||||
labels_slices[i].append(f"{E_slice_total:.1f} MeV")
|
||||
|
||||
ax.set_yscale('log')
|
||||
ax.set_ylim(0.05, 30000)
|
||||
ax.set_ylabel(slice_labels[i], fontsize=fs-6)
|
||||
ax.grid(True, which='both', ls='--', lw=0.6)
|
||||
|
||||
# Summe 0-10 cm
|
||||
ax = axes[5]
|
||||
E_sum, _ = np.histogram(df["r"], bins=r_edges, weights=df["edep"])
|
||||
sim_profile = E_sum / (N_events * A_ring)
|
||||
E_total = np.sum(E_sum) / N_events
|
||||
E_total = min(E_total, E)
|
||||
|
||||
line = ax.step(r_centers, sim_profile, where='mid', lw=2, color=col)[0]
|
||||
handles_slices[5].append(line)
|
||||
labels_slices[5].append(f"{E_total:.1f} MeV")
|
||||
|
||||
ax.set_yscale('log')
|
||||
ax.set_ylim(0.05, 30000)
|
||||
ax.set_ylabel('Sum', fontsize=fs-6)
|
||||
ax.grid(True, which='both', ls='--', lw=0.6)
|
||||
|
||||
# Füge Subplot-Legenden (nur integrierte Energie) hinzu
|
||||
for ax, h, l in zip(axes, handles_slices, labels_slices):
|
||||
ax.legend(h, l, loc='upper right', fontsize=fs-9, frameon=True)
|
||||
|
||||
# Cosmetics
|
||||
axes[-1].set_xlabel('Radius r in cm', fontsize=fs)
|
||||
for ax in axes:
|
||||
ax.set_xlim(0, r_max_cm)
|
||||
ax.tick_params(labelsize=fs-2)
|
||||
|
||||
# Gemeinsames y-Label
|
||||
fig.text(0.04, 0.5, r"$\left\langle \int \frac{dE}{dz\,dA} dz \right\rangle$ in MeV/cm$^2$",
|
||||
va='center', ha='center', rotation='vertical', fontsize=fs)
|
||||
|
||||
# Überschrift
|
||||
fig.suptitle("Transverse shower profiles in BGO (Geant4 Simulation)", fontsize=fs+2, y=0.94)
|
||||
|
||||
# Startenergie-Legende oben zwischen Titel und oberstem Plot
|
||||
fig.legend(energy_lines, labels, loc='upper center', bbox_to_anchor=(0.5,0.92),
|
||||
ncol=5, fontsize=fs-10, frameon=True)
|
||||
|
||||
plt.tight_layout(rect=[0,0,1,0.85])
|
||||
plt.savefig("plots/G4_transverse_layers_sum2.pdf", dpi=300)
|
||||
plt.show()
|
||||
|
|
@ -1,13 +1,17 @@
|
|||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
transversal shower profile theory slicewise
|
||||
transversal shower profile theory slicewise (BGO, NKG-corrected)
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from scipy.stats import gamma
|
||||
import scienceplots
|
||||
from matplotlib.lines import Line2D
|
||||
|
||||
plt.style.use(['science'])
|
||||
fs = 26
|
||||
|
||||
# --------------------------
|
||||
# Material & Parameter
|
||||
|
|
@ -20,8 +24,6 @@ Rm = 2.26
|
|||
energies = [100, 200, 500, 1000, 2000] # MeV
|
||||
colors = ["C0","C1","C2","C3","C4"]
|
||||
|
||||
fs = 24 # globale Schriftgröße
|
||||
|
||||
# --------------------------
|
||||
# Profile
|
||||
# --------------------------
|
||||
|
|
@ -30,15 +32,41 @@ def dEdz(z, E):
|
|||
a = 1 + b*(np.log(E/Ec)-0.5)
|
||||
return E * gamma.pdf(t, a, scale=1/b) / X0 # MeV/cm
|
||||
|
||||
def rho_r(r):
|
||||
return np.exp(-r**2/(2*Rm**2)) / (2*np.pi*Rm**2)
|
||||
# --------------------------
|
||||
# NKG transverse profile (safe at r=0)
|
||||
# --------------------------
|
||||
s = 1.0
|
||||
|
||||
# Verwende r_norm für Normierung (0–8 cm)
|
||||
r_norm = np.linspace(0, 2.5, 300)
|
||||
|
||||
def rho_r(r, r_norm=r_norm):
|
||||
"""
|
||||
NKG transverse profile
|
||||
normalized: ∫ 2π r ρ(r) dr = 1
|
||||
avoids divide-by-zero at r=0
|
||||
"""
|
||||
x = r / Rm
|
||||
core = np.zeros_like(r)
|
||||
mask = r > 0
|
||||
core[mask] = (x[mask]**(s-2)) * ((1 + x[mask])**(s-4.5))
|
||||
if r[0] == 0:
|
||||
core[0] = ((1e-9 / Rm)**(s-2)) * ((1 + 1e-9 / Rm)**(s-4.5))
|
||||
# Normiere über r_norm (0–8 cm)
|
||||
x_norm = r_norm / Rm
|
||||
core_norm = np.zeros_like(r_norm)
|
||||
mask_norm = r_norm > 0
|
||||
core_norm[mask_norm] = (x_norm[mask_norm]**(s-2)) * ((1 + x_norm[mask_norm])**(s-4.5))
|
||||
core_norm[0] = ((1e-9 / Rm)**(s-2)) * ((1 + 1e-9 / Rm)**(s-4.5))
|
||||
norm = np.trapz(2*np.pi*r_norm*core_norm, r_norm)
|
||||
return core / norm
|
||||
|
||||
# --------------------------
|
||||
# Gitter
|
||||
# --------------------------
|
||||
z_grid = np.linspace(0, 10, 400)
|
||||
dz = z_grid[1]-z_grid[0]
|
||||
r_grid = np.linspace(0, 8, 300)
|
||||
r_grid = np.linspace(0, 2.5, 300) # nur bis 2.5 cm
|
||||
dr = r_grid[1]-r_grid[0]
|
||||
|
||||
# --------------------------
|
||||
|
|
@ -50,24 +78,13 @@ n_layers = len(layers)
|
|||
# --------------------------
|
||||
# Figure
|
||||
# --------------------------
|
||||
plt.style.use(['science'])
|
||||
fig, axes = plt.subplots(n_layers+1, 1, figsize=(10,12), sharex=True)
|
||||
plt.subplots_adjust(hspace=0.1)
|
||||
fig.suptitle("Transverse energy deposition in BGO layers (theory)", fontsize=fs+2, y=0.95)
|
||||
|
||||
# Dummy-Linien für Startenergie-Legende
|
||||
dummy_lines = [axes[0].plot([], [], color=c, linewidth=2)[0] for c in colors]
|
||||
# Legende sichtbar unter dem Titel, innerhalb der Figure
|
||||
axes[0].legend(dummy_lines, [f"{E} MeV" for E in energies],
|
||||
ncol=5, fontsize=fs-10, frameon=True, framealpha=0.85,
|
||||
loc='upper center', bbox_to_anchor=(0.5, 1.05))
|
||||
|
||||
# Gemeinsames Y-Label
|
||||
fig.text(0.02, 0.5, r"$\langle \mathrm{d}E/(\mathrm{d}z\,\mathrm{d}A) \rangle$ in MeV/cm$^2$/nuc", va='center', rotation='vertical', fontsize=fs)
|
||||
|
||||
# --------------------------
|
||||
# Plots
|
||||
# --------------------------
|
||||
|
||||
for i, (z_min, z_max) in enumerate(layers + [(0,10)]):
|
||||
ax = axes[i]
|
||||
z_mask = (z_grid >= z_min) & (z_grid < z_max)
|
||||
|
|
@ -80,9 +97,9 @@ for i, (z_min, z_max) in enumerate(layers + [(0,10)]):
|
|||
linestyle = '--' if i==n_layers else '-' # Summe gestrichelt
|
||||
ax.plot(r_grid, radial_profile, color=color, linewidth=2, linestyle=linestyle, label=f"{E_layer:.1f} MeV")
|
||||
|
||||
ax.set_xlim(0, 8)
|
||||
ax.set_ylim(0, None)
|
||||
ax.grid(True)
|
||||
ax.set_xlim(0, 2.5)
|
||||
ax.set_ylim(0.01, 1e5)
|
||||
ax.grid(True, which='both', linestyle='--', linewidth=0.5)
|
||||
ax.tick_params(labelsize=fs-2)
|
||||
|
||||
if i < n_layers:
|
||||
|
|
@ -91,8 +108,28 @@ for i, (z_min, z_max) in enumerate(layers + [(0,10)]):
|
|||
ax.set_ylabel("Sum", fontsize=fs)
|
||||
ax.set_xlabel("Radius r in cm", fontsize=fs)
|
||||
|
||||
ax.set_yscale("log") # logarithmische y-Achse
|
||||
ax.legend(fontsize=12, frameon=True, framealpha=0.85, facecolor="white")
|
||||
|
||||
# Cosmetics
|
||||
axes[-1].set_xlabel('Radius r in cm', fontsize=fs)
|
||||
for ax in axes:
|
||||
ax.set_xlim(0,r_max_cm)
|
||||
ax.tick_params(labelsize=fs-2)
|
||||
|
||||
# Gemeinsames y-Label
|
||||
fig.text(0.04, 0.5,
|
||||
r"$\left\langle \int \frac{dE}{dz\,dA} dz \right\rangle$ [MeV/cm$^2$]",
|
||||
va='center', ha='center', rotation='vertical', fontsize=fs)
|
||||
|
||||
# Überschrift
|
||||
fig.suptitle("Transverse shower profiles in BGO (Theory)", fontsize=fs+2, y=0.98)
|
||||
|
||||
# Startenergie-Legende oberhalb oberstem Plot
|
||||
energy_lines = [Line2D([0],[0], color=c,lw=2) for c in colors]
|
||||
fig.legend(energy_lines, labels, loc='upper center', bbox_to_anchor=(0.5,0.95),
|
||||
ncol=5, fontsize=fs-10, frameon=True)
|
||||
|
||||
plt.tight_layout(rect=[0.05,0.03,0.97,0.93])
|
||||
plt.savefig("plots/shower_transverse_layers_sum.pdf")
|
||||
plt.savefig("plots/shower_transverse_layers_sum.png", dpi=300)
|
||||
|
|
|
|||
133
shower_trans_theory2.py
Normal file
133
shower_trans_theory2.py
Normal file
|
|
@ -0,0 +1,133 @@
|
|||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Transverse shower profiles in BGO
|
||||
Nur Theorie
|
||||
Legenden: obere Legende = Startenergien, rechts pro Subplot = integrierte Energie pro Slice
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from scipy.stats import gamma
|
||||
import scienceplots
|
||||
from matplotlib.lines import Line2D
|
||||
|
||||
plt.style.use(['science'])
|
||||
fs = 26
|
||||
|
||||
# Material
|
||||
X0 = 1.118
|
||||
Ec = 10.5
|
||||
b = 0.5
|
||||
Rm = 2.26
|
||||
|
||||
# Energies & Farben
|
||||
energies = [100, 200, 500, 1000, 2000]
|
||||
labels = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
|
||||
colors = ["C0", "C1", "C2", "C3", "C4"]
|
||||
|
||||
# Theorie-Funktionen
|
||||
def tmax(E):
|
||||
return np.log(E/Ec)-0.5
|
||||
|
||||
def dEdz(z, E):
|
||||
t = z/X0
|
||||
a = 1 + b*tmax(E)
|
||||
return E*gamma.pdf(t, a, scale=1/b) # MeV/cm
|
||||
|
||||
def rho_r(r):
|
||||
s = 1.0
|
||||
x = r/Rm
|
||||
core = np.zeros_like(r)
|
||||
mask = r > 0
|
||||
core[mask] = (x[mask]**(s-2)) * ((1+x[mask])**(s-4.5))
|
||||
if r[0] == 0:
|
||||
core[0] = ((1e-9/Rm)**(s-2)) * ((1+1e-9/Rm)**(s-4.5))
|
||||
norm = np.trapz(2*np.pi*r*core, r)
|
||||
return core/norm
|
||||
|
||||
# Radial & Slices
|
||||
r_max_cm = 2.5
|
||||
n_bins = 80
|
||||
r_edges = np.linspace(0, r_max_cm, n_bins+1)
|
||||
r_centers = 0.5*(r_edges[:-1] + r_edges[1:])
|
||||
r_in = r_edges[:-1]
|
||||
r_out = r_edges[1:]
|
||||
A_ring = np.pi*(r_out**2 - r_in**2)
|
||||
|
||||
z_slices = [(0,2), (2,4), (4,6), (6,8), (8,10)]
|
||||
slice_labels = ["0-2 cm", "2-4 cm", "4-6 cm", "6-8 cm", "8-10 cm"]
|
||||
|
||||
z_grid = np.linspace(0, 10, 400)
|
||||
dz = z_grid[1] - z_grid[0]
|
||||
rho = rho_r(r_centers)
|
||||
|
||||
# Figure
|
||||
fig, axes = plt.subplots(6, 1, figsize=(13,16), sharex=True,
|
||||
gridspec_kw={"height_ratios":[1,1,1,1,1,1], "hspace":0.05})
|
||||
|
||||
# Listen für Subplot-Legenden vorbereiten
|
||||
handles_slices = [[] for _ in range(len(z_slices)+1)]
|
||||
labels_slices = [[] for _ in range(len(z_slices)+1)]
|
||||
|
||||
for E, lab, col in zip(energies, labels, colors):
|
||||
for i, (z1, z2) in enumerate(z_slices):
|
||||
ax = axes[i]
|
||||
mask = (z_grid >= z1) & (z_grid < z2)
|
||||
theory_profile = np.sum(dEdz(z_grid[mask], E)[:,None]*rho[None,:]*dz, axis=0)
|
||||
|
||||
# Slice-Energie
|
||||
E_slice_total = np.sum(theory_profile*A_ring)
|
||||
E_slice_total = min(E_slice_total, E)
|
||||
|
||||
line, = ax.plot(r_centers, theory_profile, lw=2, color=col)
|
||||
|
||||
handles_slices[i].append(line)
|
||||
labels_slices[i].append(f"{E_slice_total:.1f} MeV") # nur integrierte Energie
|
||||
|
||||
ax.set_yscale('log')
|
||||
ax.set_ylim(0.05, 30000)
|
||||
ax.set_ylabel(slice_labels[i], fontsize=fs-6)
|
||||
ax.grid(True, which='both', ls='--', lw=0.6)
|
||||
|
||||
# Summe 0-10 cm
|
||||
ax = axes[5]
|
||||
# Theorie Slice-wise
|
||||
mask = (z_grid >= z1) & (z_grid < z2)
|
||||
theory_profile = np.sum(dEdz(z_grid[mask], E)[:, None] * rho[None, :] * dz, axis=0)
|
||||
E_total = np.sum(theory_profile*A_ring)
|
||||
E_total = min(E_total, E)
|
||||
|
||||
line, = ax.plot(r_centers, theory_profile, lw=2, color=col)
|
||||
handles_slices[5].append(line)
|
||||
labels_slices[5].append(f"{E_total:.1f} MeV")
|
||||
|
||||
ax.set_yscale('log')
|
||||
ax.set_ylim(0.05, 30000)
|
||||
ax.set_ylabel('Sum', fontsize=fs-6)
|
||||
ax.grid(True, which='both', ls='--', lw=0.6)
|
||||
|
||||
# Subplot-Legenden (integrierte Energie pro Slice)
|
||||
for ax, h, l in zip(axes, handles_slices, labels_slices):
|
||||
ax.legend(h, l, loc='upper right', fontsize=fs-9, frameon=True)
|
||||
|
||||
# Cosmetics
|
||||
axes[-1].set_xlabel('Radius r in cm', fontsize=fs)
|
||||
for ax in axes:
|
||||
ax.set_xlim(0, r_max_cm)
|
||||
ax.tick_params(labelsize=fs-2)
|
||||
|
||||
fig.text(0.04, 0.5,
|
||||
r"$\left\langle \int \frac{dE}{dz\,dA} dz \right\rangle$ in MeV/cm$^2$",
|
||||
va='center', ha='center', rotation='vertical', fontsize=fs)
|
||||
|
||||
fig.suptitle("Transverse shower profiles in BGO (Theory)", fontsize=fs+2, y=0.94)
|
||||
|
||||
# Startenergie-Legende oben zwischen Titel und oberstem Plot
|
||||
energy_lines = [Line2D([0],[0], color=c, lw=2) for c in colors]
|
||||
fig.legend(energy_lines, labels, loc='upper center', bbox_to_anchor=(0.5,0.92),
|
||||
ncol=5, fontsize=fs-10, frameon=True)
|
||||
|
||||
plt.tight_layout(rect=[0,0,1,0.85])
|
||||
plt.savefig("plots/shower_transverse_layers_sum2.pdf", dpi=300)
|
||||
plt.show()
|
||||
222
showermap_comparison.py
Normal file
222
showermap_comparison.py
Normal file
|
|
@ -0,0 +1,222 @@
|
|||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
|
||||
"""
|
||||
6-Panel plot:
|
||||
- Top-left: Theory 2D contours (absolute dE/dz in MeV/cm²)
|
||||
- Other 5: Geant4 heatmaps for 100,200,500,1GeV,2GeV
|
||||
- Shared colorbar for Geant4 heatmaps, vmax=0.8 MeV/cm²
|
||||
- Subcaptions: "Theory", "100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"
|
||||
- Theory contours overlayed on Geant4 heatmaps
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib.gridspec import GridSpec
|
||||
from matplotlib.lines import Line2D
|
||||
from scipy.stats import gamma
|
||||
import scienceplots
|
||||
|
||||
plt.style.use(['science'])
|
||||
fs = 22
|
||||
|
||||
# --------------------------
|
||||
# Material & BGO parameters
|
||||
# --------------------------
|
||||
X0_cm = 1.118
|
||||
Ec_MeV = 10.5
|
||||
b = 0.5
|
||||
Rm_cm = 2.26
|
||||
max_depth_cm = 10
|
||||
|
||||
# --------------------------
|
||||
# Energies
|
||||
# --------------------------
|
||||
energies_MeV = [100, 200, 500, 1000, 2000]
|
||||
colors = ["C0", "C1", "C2", "C3", "C4"]
|
||||
|
||||
# --------------------------
|
||||
# Theory functions
|
||||
# --------------------------
|
||||
def tmax(E, Ec):
|
||||
return np.log(E / Ec) - 0.5
|
||||
|
||||
def dEdz(z, E, X0=X0_cm, b=0.5, Ec=Ec_MeV):
|
||||
"""Longitudinal deposition MeV/cm"""
|
||||
t = z / X0
|
||||
a = 1 + b * tmax(E, Ec)
|
||||
return E * gamma.pdf(t, a, scale=1 / b) # MeV/cm
|
||||
|
||||
s = 1.0
|
||||
def rho_r(r, Rm=Rm_cm, s_val=1.0):
|
||||
x = r / Rm
|
||||
core = np.zeros_like(r)
|
||||
mask = r > 0
|
||||
core[mask] = (x[mask] ** (s_val - 2)) * ((1 + x[mask]) ** (s_val - 4.5))
|
||||
if r[0] == 0:
|
||||
core[0] = ((1e-9 / Rm) ** (s_val - 2)) * ((1 + 1e-9 / Rm) ** (s_val - 4.5))
|
||||
norm = np.trapz(2 * np.pi * r * core, r)
|
||||
return core / norm
|
||||
|
||||
# --------------------------
|
||||
# Radial & depth grid
|
||||
# --------------------------
|
||||
r_max = 2.5
|
||||
r = np.linspace(0, r_max, 300)
|
||||
z = np.linspace(0, max_depth_cm, 400)
|
||||
Z, R = np.meshgrid(z, r, indexing="ij")
|
||||
|
||||
# --------------------------
|
||||
# Geant4 files
|
||||
# --------------------------
|
||||
files = [
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_100MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_200MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_500MeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_1GeV_0.hits",
|
||||
"/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_2GeV_0.hits",
|
||||
]
|
||||
|
||||
# --------------------------
|
||||
# Figure layout
|
||||
# --------------------------
|
||||
fig = plt.figure(figsize=(18, 10))
|
||||
gs = GridSpec(2, 3, figure=fig, wspace=0.3, hspace=0.4)
|
||||
|
||||
axes = [
|
||||
fig.add_subplot(gs[0, 0]), # Theory contours
|
||||
fig.add_subplot(gs[0, 1]), # 100 MeV
|
||||
fig.add_subplot(gs[0, 2]), # 200 MeV
|
||||
fig.add_subplot(gs[1, 0]), # 500 MeV
|
||||
fig.add_subplot(gs[1, 1]), # 1 GeV
|
||||
fig.add_subplot(gs[1, 2]), # 2 GeV
|
||||
]
|
||||
|
||||
subtitles = ["Theory", "100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
|
||||
|
||||
# --------------------------
|
||||
# 1) Theory contour plot
|
||||
# --------------------------
|
||||
ax = axes[0]
|
||||
heatmap_levels = [0.1, 1, 10] # MeV/cm²
|
||||
linestyles = [":", "--", "-"]
|
||||
|
||||
for E, col in zip(energies_MeV, colors):
|
||||
long_prof = dEdz(z, E)
|
||||
shower_2D = np.zeros((len(z), len(r)))
|
||||
rho = rho_r(r)
|
||||
rho[0] = rho[1]
|
||||
for i in range(len(z)):
|
||||
shower_2D[i, :] = long_prof[i] * rho
|
||||
|
||||
ax.contour(
|
||||
R, Z, shower_2D,
|
||||
levels=heatmap_levels,
|
||||
colors=[col] * len(heatmap_levels),
|
||||
linestyles=linestyles,
|
||||
linewidths=1.5
|
||||
)
|
||||
|
||||
ax.set_xlabel("Radius r in cm", fontsize=fs)
|
||||
ax.set_ylabel("Depth z in cm", fontsize=fs)
|
||||
ax.set_xlim(0, r_max)
|
||||
ax.set_ylim(max_depth_cm, 0)
|
||||
ax.set_title("2D Energy density contours", fontsize=fs)
|
||||
ax.grid(True, linestyle="--", linewidth=0.5)
|
||||
ax.tick_params(labelsize=fs)
|
||||
|
||||
line_legend = [Line2D([0], [0], color="k", ls=ls, lw=2) for ls in linestyles]
|
||||
ax.legend(
|
||||
line_legend,
|
||||
[f"{lvl} MeV/cm²" for lvl in heatmap_levels],
|
||||
title="Contour level",
|
||||
fontsize=fs - 2,
|
||||
title_fontsize=fs,
|
||||
loc="lower right",
|
||||
frameon=True,
|
||||
)
|
||||
|
||||
# --------------------------
|
||||
# 2–6) Geant4 heatmaps + Theory overlay
|
||||
# --------------------------
|
||||
vmin_cb = 0.0
|
||||
vmax_cb = 0.8
|
||||
im_list = []
|
||||
|
||||
for ax, file, subtitle, col in zip(axes[1:], files, subtitles[1:], colors):
|
||||
df = pd.read_csv(file, sep="\t")
|
||||
df = df[df["r"] <= r_max]
|
||||
|
||||
dz = 0.1
|
||||
dr = 0.1
|
||||
z_edges = np.arange(0, max_depth_cm + dz, dz)
|
||||
r_edges = np.arange(0, r_max + dr, dr)
|
||||
|
||||
H, _, _ = np.histogram2d(
|
||||
df["z"], df["r"],
|
||||
bins=[z_edges, r_edges],
|
||||
weights=df["edep"]
|
||||
)
|
||||
counts, _, _ = np.histogram2d(
|
||||
df["z"], df["r"],
|
||||
bins=[z_edges, r_edges]
|
||||
)
|
||||
|
||||
H_avg = np.zeros_like(H)
|
||||
mask = counts > 0
|
||||
H_avg[mask] = H[mask] / counts[mask]
|
||||
|
||||
Z_grid, R_grid = np.meshgrid(
|
||||
(z_edges[:-1] + z_edges[1:]) / 2,
|
||||
(r_edges[:-1] + r_edges[1:]) / 2,
|
||||
indexing="ij"
|
||||
)
|
||||
|
||||
im = ax.pcolormesh(
|
||||
R_grid, Z_grid, H_avg,
|
||||
shading="auto",
|
||||
cmap="viridis",
|
||||
vmin=vmin_cb,
|
||||
vmax=vmax_cb
|
||||
)
|
||||
im_list.append(im)
|
||||
|
||||
# Theory overlay contours
|
||||
E_val = float(subtitle.split()[0])
|
||||
long_prof = dEdz(z, E_val)
|
||||
shower_2D = np.zeros((len(z), len(r)))
|
||||
rho = rho_r(r)
|
||||
rho[0] = rho[1]
|
||||
for i in range(len(z)):
|
||||
shower_2D[i, :] = long_prof[i] * rho
|
||||
|
||||
ax.contour(
|
||||
R, Z, shower_2D,
|
||||
levels=heatmap_levels,
|
||||
colors="white",
|
||||
linestyles=linestyles,
|
||||
linewidths=1
|
||||
)
|
||||
|
||||
ax.set_title(subtitle, fontsize=fs, color=col, pad=10)
|
||||
ax.set_xlabel("Radius r in cm", fontsize=fs)
|
||||
ax.set_ylabel("Depth z in cm", fontsize=fs)
|
||||
ax.set_ylim(max_depth_cm, 0)
|
||||
ax.set_xlim(0, r_max)
|
||||
ax.grid(True, linestyle="--", linewidth=0.5)
|
||||
ax.tick_params(labelsize=fs)
|
||||
|
||||
# --------------------------
|
||||
# Shared colorbar
|
||||
# --------------------------
|
||||
ymax = axes[1].get_position().y1
|
||||
ymin = axes[-1].get_position().y0
|
||||
cbar_ax = fig.add_axes([0.93, ymin, 0.015, ymax - ymin])
|
||||
cbar = fig.colorbar(im_list[0], cax=cbar_ax)
|
||||
cbar.set_label("dE/dz in MeV/cm²", fontsize=fs)
|
||||
cbar.ax.tick_params(labelsize=fs - 2)
|
||||
|
||||
plt.tight_layout(rect=[0, 0, 0.92, 0.95])
|
||||
plt.savefig("plots/shower_map_comparison.pdf", dpi=300)
|
||||
plt.show()
|
||||
Loading…
Add table
Add a link
Reference in a new issue