Compare commits
2 commits
7e703cd499
...
233684286d
| Author | SHA1 | Date | |
|---|---|---|---|
|
233684286d | ||
|
|
8affb23be8 |
12 changed files with 262 additions and 191 deletions
|
|
@ -85,3 +85,68 @@ def photon_emission_ava(n,E,m_0,Z,lam1,lam2):
|
|||
#else: print('no: E=',E[i],";v=",v, "n=")
|
||||
|
||||
return result
|
||||
|
||||
|
||||
def cherenkov_theta(n, E_kin, m):
|
||||
"""
|
||||
Calculate Cherenkov angle theta_C as a function of kinetic energy.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
n : float
|
||||
Refractive index
|
||||
E_kin : array_like
|
||||
Kinetic energy [J]
|
||||
m : float
|
||||
Particle mass [kg]
|
||||
|
||||
Returns
|
||||
-------
|
||||
theta : ndarray
|
||||
Cherenkov angle [rad]; zero below threshold
|
||||
"""
|
||||
|
||||
gamma = 1.0 + E_kin / (m * c**2)
|
||||
beta2 = 1.0 - 1.0 / gamma**2
|
||||
beta2 = np.clip(beta2, 0.0, None)
|
||||
beta = np.sqrt(beta2)
|
||||
|
||||
cos_theta = 1.0 / (beta * n)
|
||||
|
||||
theta = np.zeros_like(E_kin)
|
||||
mask = cos_theta < 1.0 # Cherenkov condition
|
||||
theta[mask] = np.arccos(cos_theta[mask])
|
||||
|
||||
return theta
|
||||
|
||||
def cherenkov_theta_deg(n, E_kin, m):
|
||||
"""
|
||||
Cherenkov angle in degrees as a function of kinetic energy.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
n : float
|
||||
Refractive index
|
||||
E_kin : array_like
|
||||
Kinetic energy [J]
|
||||
m : float
|
||||
Particle mass [kg]
|
||||
|
||||
Returns
|
||||
-------
|
||||
theta_deg : ndarray
|
||||
Cherenkov angle [deg]; zero below threshold
|
||||
"""
|
||||
|
||||
gamma = 1.0 + E_kin / (m * c**2)
|
||||
beta2 = 1.0 - 1.0 / gamma**2
|
||||
beta2 = np.clip(beta2, 0.0, None)
|
||||
beta = np.sqrt(beta2)
|
||||
|
||||
cos_theta = 1.0 / (beta * n)
|
||||
|
||||
theta = np.zeros_like(E_kin)
|
||||
mask = cos_theta < 1.0
|
||||
theta[mask] = np.arccos(cos_theta[mask])
|
||||
|
||||
return np.degrees(theta)
|
||||
|
|
@ -65,6 +65,7 @@ n_p = photon_emission_ava(n,E0,m_p,z_p,lam1,lam2)
|
|||
|
||||
|
||||
fs = 24
|
||||
fs = 32
|
||||
|
||||
#plt.text(1, 76.5, '76.5', color='black', fontsize=fs-4, verticalalignment='bottom', horizontalalignment='left')
|
||||
#plt.text(1, 220, '306.1', color='green', fontsize=fs-4, verticalalignment='bottom', horizontalalignment='left')
|
||||
|
|
@ -109,6 +110,65 @@ plt.legend(fontsize=fs-1, loc="lower right",ncol=1,
|
|||
plt.savefig("plots/frank-tamm")
|
||||
#plt.show()
|
||||
|
||||
|
||||
##################### ANGLE #########################
|
||||
plt.figure(figsize=(15, 10))
|
||||
plt.style.use(['science'])
|
||||
|
||||
# Cherenkov angles
|
||||
theta_e = cherenkov_theta_deg(1.05, E0, m_e)
|
||||
theta_m = cherenkov_theta_deg(1.05, E0, m_mu)
|
||||
theta_p = cherenkov_theta_deg(1.05, E0, m_p)
|
||||
theta_a = cherenkov_theta_deg(1.05, E0, m_a)
|
||||
|
||||
theta_e2 = cherenkov_theta_deg(1.07, E0, m_e)
|
||||
theta_m2 = cherenkov_theta_deg(1.07, E0, m_mu)
|
||||
theta_p2 = cherenkov_theta_deg(1.07, E0, m_p)
|
||||
theta_a2 = cherenkov_theta_deg(1.07, E0, m_a)
|
||||
|
||||
# Grenzwinkel
|
||||
theta_max_105 = np.degrees(np.arccos(1 / 1.05))
|
||||
theta_max_107 = np.degrees(np.arccos(1 / 1.07))
|
||||
|
||||
# Farben exakt wie Frank–Tamm-Plot
|
||||
plt.plot(E0*J_to_MeV, theta_e, color='blue', lw=2, label='electron')
|
||||
plt.plot(E0*J_to_MeV, theta_m, color=(0.5, 0.3, 0.8), lw=2, label='muon')
|
||||
plt.plot(E0*J_to_MeV, theta_p, color='red', lw=2, label='proton')
|
||||
plt.plot(E0*J_to_MeV, theta_a, color='green', lw=2, label='helium')
|
||||
|
||||
plt.plot(E0*J_to_MeV, theta_e2, color='blue', lw=2, ls='--')
|
||||
plt.plot(E0*J_to_MeV, theta_m2, color=(0.5, 0.3, 0.8), lw=2, ls='--')
|
||||
plt.plot(E0*J_to_MeV, theta_p2, color='red', lw=2, ls='--')
|
||||
plt.plot(E0*J_to_MeV, theta_a2, color='green', lw=2, ls='--')
|
||||
|
||||
# Grenzwinkel-Linien
|
||||
plt.axhline(theta_max_105, color='black', lw=2)
|
||||
plt.axhline(theta_max_107, color='black', lw=2, ls='--')
|
||||
|
||||
plt.text(1, theta_max_105, rf'n=1.05: ${theta_max_105:.1f}^\circ$', color='black', fontsize=fs-4, verticalalignment='bottom', horizontalalignment='left')
|
||||
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.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.legend(fontsize=fs-2, loc='lower right',
|
||||
frameon=True, framealpha=0.4, fancybox=True)
|
||||
|
||||
plt.savefig("plots/cherenkov_angle_deg")
|
||||
|
||||
|
||||
###############################################################
|
||||
|
||||
N = []
|
||||
# for E in E0_p:
|
||||
# n = photon_emission_carlotta(n,E*1e-6,m_p,z_p,lam1,lam2)
|
||||
|
|
@ -194,7 +254,7 @@ fig, ax1 = plt.subplots(figsize=(15, 10))
|
|||
plt.style.use(['science'])
|
||||
|
||||
# --- linke y-Achse: Cherenkov-Spektrum ---
|
||||
ax1.plot(
|
||||
ax1.step(
|
||||
lam_centers,
|
||||
N_lambda,
|
||||
lw=2,
|
||||
|
|
@ -434,7 +494,7 @@ plt.minorticks_on()
|
|||
ax2.set_yscale('log')
|
||||
#ax1.set_yscale('log')
|
||||
|
||||
plt.savefig("plots/prediction.pdf")
|
||||
#plt.savefig("plots/prediction.pdf")
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -60,7 +60,7 @@ 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}\;[\mathrm{cm}^{-1}]$',
|
||||
r'$\frac{dN_{\mathrm{photons}}}{dx}$ in $\mathrm{cm}^{-1}$',
|
||||
fontsize=fs+2
|
||||
)
|
||||
|
||||
|
|
@ -87,6 +87,14 @@ lam_centers = Nlam[:-1] + 5
|
|||
|
||||
n_list = [1.03, 1.05, 1.07, 1.09]
|
||||
|
||||
# ---------- Farben konsistent ----------
|
||||
colors = {
|
||||
1.03: 'tab:orange',
|
||||
1.05: 'tab:blue',
|
||||
1.07: 'tab:green',
|
||||
1.09: 'tab:red'
|
||||
}
|
||||
|
||||
plt.figure(figsize=(15, 10))
|
||||
plt.style.use(['science'])
|
||||
|
||||
|
|
@ -106,17 +114,18 @@ for n_val in n_list:
|
|||
|
||||
N_lambda = np.array(N_lambda)
|
||||
|
||||
plt.plot(
|
||||
plt.step(
|
||||
lam_centers,
|
||||
N_lambda,
|
||||
lw=2,
|
||||
color=colors[n_val],
|
||||
label=fr'$n={n_val:.2f}$'
|
||||
)
|
||||
|
||||
# ---------- Plot-Format ----------
|
||||
plt.xlabel(r'Wavelength $\lambda$ [nm]', fontsize=fs+2)
|
||||
plt.xlabel(r'Wavelength $\lambda$ in nm', fontsize=fs+2)
|
||||
plt.ylabel(
|
||||
r'$\frac{dN_{\mathrm{photons}}}{dx}\;[\mathrm{cm}^{-1}]$',
|
||||
r'$\frac{dN_{\mathrm{photons}}}{dx}$ in $\mathrm{cm}^{-1}$',
|
||||
fontsize=fs+2
|
||||
)
|
||||
|
||||
|
|
@ -165,14 +174,6 @@ T_centers = T_percent_centers / 100.0
|
|||
fig, ax1 = plt.subplots(figsize=(15, 10))
|
||||
plt.style.use(['science'])
|
||||
|
||||
# ---------- Farben konsistent ----------
|
||||
colors = {
|
||||
1.03: 'tab:blue',
|
||||
1.05: 'tab:orange',
|
||||
1.07: 'tab:green',
|
||||
1.09: 'tab:red'
|
||||
}
|
||||
|
||||
# ---------- Cherenkov-Spektren ----------
|
||||
spectra = {}
|
||||
|
||||
|
|
@ -211,9 +212,9 @@ for n_val in n_list:
|
|||
)
|
||||
|
||||
# ---------- Linke Achse ----------
|
||||
ax1.set_xlabel(r'Wavelength $\lambda$ [nm]', fontsize=fs+2)
|
||||
ax1.set_xlabel(r'Wavelength $\lambda$ in nm', fontsize=fs+2)
|
||||
ax1.set_ylabel(
|
||||
r'$\frac{\mathrm{d}N_{\mathrm{photons}}}{\mathrm{d}x}\;[\mathrm{cm}^{-1}]$',
|
||||
r'$\frac{\mathrm{d}N_{\mathrm{photons}}}{\mathrm{d}x}$ in $\mathrm{cm}^{-1}$',
|
||||
fontsize=fs+2
|
||||
)
|
||||
ax1.tick_params(axis='both', labelsize=fs)
|
||||
|
|
@ -235,7 +236,7 @@ ax2.plot(
|
|||
label='Transmittance'
|
||||
)
|
||||
|
||||
ax2.set_ylabel(r'Transmittance (\%)', fontsize=fs+2)
|
||||
ax2.set_ylabel(r'Transmittance in (\%)', fontsize=fs+2)
|
||||
ax2.tick_params(axis='y', labelsize=fs)
|
||||
|
||||
# ---------- Titel ----------
|
||||
|
|
@ -324,7 +325,12 @@ 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'\;[\mathrm{cm}^{-1}]$', 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)
|
||||
|
|
@ -358,7 +364,7 @@ plt.plot(n_values, E_thr_mu, lw=2, color=(0.5, 0.3, 0.8), label='muon')
|
|||
plt.yscale('log')
|
||||
|
||||
plt.xlabel(r'Refractive index $n$', fontsize=fs+2)
|
||||
plt.ylabel(r'Cherenkov threshold energy $E_{\mathrm{thr}}$ [MeV]', fontsize=fs+2)
|
||||
plt.ylabel(r'$E_{\mathrm{thr}}$ in MeV', fontsize=fs+2)
|
||||
|
||||
plt.title(
|
||||
r'Cherenkov threshold energy',
|
||||
|
|
@ -403,8 +409,8 @@ plt.plot(
|
|||
|
||||
plt.xlabel(r'Refractive index $n$', fontsize=fs+2)
|
||||
plt.ylabel(
|
||||
r'Cherenkov threshold energy $E_{\mathrm{thr}}$ [GeV]',
|
||||
fontsize=26
|
||||
r'$E_{\mathrm{thr}}$ in GeV',
|
||||
fontsize=fs+2
|
||||
)
|
||||
|
||||
plt.title(
|
||||
|
|
|
|||
|
|
@ -9,12 +9,16 @@ import argparse
|
|||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
import scienceplots
|
||||
|
||||
plt.style.use(['science'])
|
||||
fs = 24
|
||||
|
||||
def main():
|
||||
parser = argparse.ArgumentParser(description="2D shower heatmap (dE/dz per step)")
|
||||
parser.add_argument("file", help="Input .hits or .csv file")
|
||||
parser.add_argument("--rmax", type=float, default=8.0, help="Max radius [cm]")
|
||||
parser.add_argument("--zmax", type=float, default=10.0, help="Max depth [cm]")
|
||||
parser.add_argument("--rmax", type=float, default=8.0, help="Max radius in cm")
|
||||
parser.add_argument("--zmax", type=float, default=10.0, help="Max depth in cm")
|
||||
parser.add_argument("--out", default="shower_heatmap.png", help="Output image")
|
||||
args = parser.parse_args()
|
||||
|
||||
|
|
@ -54,14 +58,16 @@ def main():
|
|||
)
|
||||
|
||||
plt.figure(figsize=(8, 6))
|
||||
im = plt.pcolormesh(R, Z, H_avg, shading="auto", cmap="inferno")
|
||||
plt.colorbar(im, label="dE/dz [MeV/cm]")
|
||||
im = plt.pcolormesh(R, Z, H_avg, shading="auto", cmap="viridis")
|
||||
cbar = plt.colorbar(im) # Colorbar erzeugen
|
||||
cbar.set_label(r"$\mathrm{d}E/\mathrm{d}z$ in MeV/cm", fontsize=fs)
|
||||
cbar.ax.tick_params(labelsize=fs-2)
|
||||
|
||||
plt.xlabel("Radius r [cm]")
|
||||
plt.ylabel("Depth z [cm]")
|
||||
plt.xlabel("Radius r in cm", fontsize=fs)
|
||||
plt.ylabel("Depth z in cm", fontsize=fs)
|
||||
plt.ylim(zmax, 0)
|
||||
plt.tick_params(axis='both', which='major', labelsize=fs-5)
|
||||
plt.title("Shower energy deposition (simulation)")
|
||||
plt.tick_params(axis='both', which='major', labelsize=fs-2)
|
||||
plt.title("Shower energy deposition 2GeV (Geant4)", fontsize=fs)
|
||||
|
||||
plt.tight_layout()
|
||||
plt.savefig(args.out, dpi=300)
|
||||
|
|
|
|||
|
|
@ -9,6 +9,7 @@ Comparable to analytic shower theory
|
|||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
import scienceplots
|
||||
|
||||
# -------------------------------------------------
|
||||
# Input files
|
||||
|
|
@ -29,12 +30,13 @@ colors = ["C0", "C1", "C2", "C3", "C4"]
|
|||
# -------------------------------------------------
|
||||
z_max_cm = 10.0
|
||||
n_bins = 100
|
||||
fs = 18
|
||||
fs = 26
|
||||
|
||||
# -------------------------------------------------
|
||||
# Plot
|
||||
# -------------------------------------------------
|
||||
plt.figure(figsize=(12, 6))
|
||||
plt.style.use(['science'])
|
||||
|
||||
for file, label, color in zip(files, labels, colors):
|
||||
|
||||
|
|
@ -63,74 +65,14 @@ for file, label, color in zip(files, labels, colors):
|
|||
label=label
|
||||
)
|
||||
|
||||
plt.xlabel("Depth z [cm]", fontsize=fs)
|
||||
plt.ylabel(r"$\langle dE/dz \rangle$ [MeV/cm/nuc]", fontsize=fs)
|
||||
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.tick_params(labelsize=fs-2)
|
||||
plt.legend(title="Initial Energy")
|
||||
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()
|
||||
|
||||
|
||||
# # --------- 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 ---------
|
||||
# z_max_cm = 10.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")
|
||||
|
||||
# z_edges = np.linspace(0, z_max_cm, n_bins + 1)
|
||||
# dz = z_edges[1] - z_edges[0]
|
||||
# z_centers = (z_edges[:-1] + z_edges[1:]) / 2
|
||||
|
||||
# # Energie pro z-Bin
|
||||
# E_sum, _ = np.histogram(df["z"], bins=z_edges, weights=df["edep"])
|
||||
# counts, _ = np.histogram(df["z"], bins=z_edges)
|
||||
|
||||
# # Mittelwert pro Step → dE/dz [MeV/cm]
|
||||
# profile = np.zeros_like(E_sum)
|
||||
# mask = counts > 0
|
||||
# profile[mask] = E_sum[mask] / counts[mask] / dz
|
||||
|
||||
# # Step-Plot (physikalisch korrekt für Histogramme)
|
||||
# plt.step(
|
||||
# z_centers,
|
||||
# profile,
|
||||
# where="mid",
|
||||
# lw=2,
|
||||
# color=color,
|
||||
# label=label
|
||||
# )
|
||||
|
||||
# plt.xlabel("Depth z [cm]",fontsize=fs)
|
||||
# plt.ylabel("dE/dz [MeV/cm]",fontsize=fs)
|
||||
# plt.title("Longitudinal shower profile", fontsize=fs+2)
|
||||
# plt.xlim(0, z_max_cm)
|
||||
# plt.tick_params(axis='both', which='major', labelsize=fs-3)
|
||||
# plt.grid(True, ls="--", lw=0.5)
|
||||
# plt.legend(title="Initial Energy")
|
||||
|
||||
# plt.tight_layout()
|
||||
# plt.savefig("plots/shower_longitudinal.png", dpi=300)
|
||||
# plt.show()
|
||||
|
|
|
|||
|
|
@ -1,6 +1,13 @@
|
|||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Longitudinal energy deposition in BGO by radial ring (theory)
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from scipy.stats import gamma
|
||||
import scienceplots
|
||||
|
||||
# --------------------------
|
||||
# Material & Parameter
|
||||
|
|
@ -13,6 +20,8 @@ Rm = 2.26
|
|||
energies = [100, 200, 500, 1000, 2000] # MeV
|
||||
colors = ["C0","C1","C2","C3","C4"]
|
||||
|
||||
fs = 24 # globale Schriftgröße
|
||||
|
||||
# --------------------------
|
||||
# Profile
|
||||
# --------------------------
|
||||
|
|
@ -41,18 +50,13 @@ n_rings = len(rings)
|
|||
# --------------------------
|
||||
# Figure
|
||||
# --------------------------
|
||||
plt.style.use(['science'])
|
||||
fig, axes = plt.subplots(n_rings+1, 1, figsize=(10,12), sharex=True)
|
||||
plt.subplots_adjust(hspace=0.1)
|
||||
fig.suptitle("Longitudinal energy deposition in BGO by radial ring (theory)", fontsize=16, y=0.95)
|
||||
|
||||
# Dummy-Linien für Startenergie-Legende (nur eine Zeile, 5 Spalten)
|
||||
# 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=12, frameon=True, framealpha=0.85, facecolor="white",
|
||||
# loc='upper center', bbox_to_anchor=(0.5, 1.12))
|
||||
fig.suptitle("Longitudinal energy deposition in BGO by radial ring (theory)", fontsize=fs+2, y=0.95)
|
||||
|
||||
# Gemeinsames Y-Label
|
||||
fig.text(0.02, 0.5, r"$\langle dE/dz \rangle$ [MeV/cm]", va='center', rotation='vertical', fontsize=16)
|
||||
fig.text(0.02, 0.5, r"$\langle \mathrm{d}E/\mathrm{d}z \rangle$ in MeV/cm", va='center', rotation='vertical', fontsize=fs)
|
||||
|
||||
# --------------------------
|
||||
# Plots
|
||||
|
|
@ -64,7 +68,9 @@ for i, (r_min, r_max) in enumerate(rings + [(0,8)]):
|
|||
for E, color in zip(energies, colors):
|
||||
long_prof = dEdz(z_grid, E)
|
||||
radial_profile_2d = long_prof[:, None] * rho_r(r_grid)[None, :]
|
||||
# Integration über Ringfläche
|
||||
ring_profile = np.trapz(radial_profile_2d[:, r_mask]*2*np.pi*r_grid[r_mask], x=r_grid[r_mask], axis=1)
|
||||
# Gesamtenergie im Ring
|
||||
E_ring = np.trapz(ring_profile, x=z_grid)
|
||||
linestyle = '--' if i==n_rings else '-'
|
||||
ax.plot(z_grid, ring_profile, color=color, linewidth=2, linestyle=linestyle, label=f"{E_ring:.1f} MeV")
|
||||
|
|
@ -72,15 +78,15 @@ for i, (r_min, r_max) in enumerate(rings + [(0,8)]):
|
|||
ax.set_xlim(0, 10)
|
||||
ax.set_ylim(0, None)
|
||||
ax.grid(True)
|
||||
ax.tick_params(labelsize=14)
|
||||
ax.tick_params(labelsize=fs-2)
|
||||
|
||||
if i < n_rings:
|
||||
ax.set_ylabel(f"r={r_min}-{r_max} cm", fontsize=14)
|
||||
ax.set_ylabel(f"r={r_min}-{r_max} cm", fontsize=fs)
|
||||
else:
|
||||
ax.set_ylabel("Sum", fontsize=14)
|
||||
ax.set_xlabel("Depth z [cm]", fontsize=14)
|
||||
ax.set_ylabel("Sum", fontsize=fs)
|
||||
ax.set_xlabel("Depth z in cm", fontsize=fs)
|
||||
|
||||
ax.legend(fontsize=12, frameon=True, framealpha=0.85, facecolor="white", loc='upper right')
|
||||
ax.legend(fontsize=fs-10, frameon=True, framealpha=0.85, facecolor="white", loc='upper right')
|
||||
|
||||
plt.tight_layout(rect=[0.05,0.03,0.97,0.93])
|
||||
plt.savefig("plots/shower_longitudinal_rings_sum.pdf")
|
||||
|
|
|
|||
|
|
@ -1,14 +1,15 @@
|
|||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
overwiev: 2D map, transversal and longitudinal shower profile theory
|
||||
Overview: 2D map, transversal and longitudinal shower profile theory
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib.gridspec import GridSpec
|
||||
from scipy.stats import gamma, norm
|
||||
from scipy.stats import gamma
|
||||
from matplotlib.lines import Line2D
|
||||
import scienceplots
|
||||
|
||||
# --------------------------
|
||||
# Funktionen
|
||||
|
|
@ -42,7 +43,7 @@ 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 = 16 # Schriftgröße
|
||||
fs = 22 # Schriftgröße
|
||||
|
||||
# Heatmap Konturen in Prozent
|
||||
heatmap_levels = [1, 10, 50, 90]
|
||||
|
|
@ -59,6 +60,7 @@ Z, R = np.meshgrid(z, r, indexing="ij")
|
|||
# Figure und GridSpec
|
||||
# --------------------------
|
||||
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])
|
||||
|
|
@ -73,14 +75,13 @@ for E, color in zip(energies_MeV, colors):
|
|||
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))
|
||||
#trans_prof = np.exp(-(R**2)/(2*sigma_r**2)) / (2*np.pi*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])
|
||||
ax_heat.contour(R, Z, shower_norm, levels=[pct], colors=[color], linestyles=[ls], linewidths=2)
|
||||
|
||||
ax_heat.grid(True, linestyle='--', linewidth=0.5)
|
||||
ax_heat.set_ylabel("Depth z [cm]", fontsize=fs)
|
||||
ax_heat.set_ylabel("Depth z in cm", fontsize=fs)
|
||||
ax_heat.set_ylim(max_depth_cm, 0)
|
||||
ax_heat.set_xlim(0,8)
|
||||
ax_heat.set_title("2D contourlines", fontsize=fs)
|
||||
|
|
@ -88,7 +89,7 @@ ax_heat.tick_params(labelsize=fs)
|
|||
|
||||
# Molière Linien (Heatmap)
|
||||
molier_radii = [Rm_cm, 2*Rm_cm]
|
||||
molier_widths = [1.5, 2.0]
|
||||
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)
|
||||
|
|
@ -102,63 +103,35 @@ for E, color in zip(energies_MeV, colors):
|
|||
ax_long.plot(profile, t*X0_cm, color=color, linewidth=2, label=f"{E/1000:.1f} GeV")
|
||||
|
||||
ax_long.set_title("Longitudinal Profiles", fontsize=fs)
|
||||
ax_long.set_xlabel("dE/dz [MeV/cm]", 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_ylim(max_depth_cm, 0)
|
||||
ax_long.tick_params(labelsize=fs)
|
||||
#ax_long.legend(title="Initial Energy", fontsize=fs-2, title_fontsize=fs, loc="upper right", frameon=True)
|
||||
|
||||
# --------------------------
|
||||
# Transversale Profile (korrekt normiert wie Plot 1.1)
|
||||
# Transversale Profile (integrated over depth)
|
||||
# --------------------------
|
||||
energy_lines = []
|
||||
|
||||
#lokal im Schauermaximum
|
||||
# for E, color in zip(energies_MeV, colors):
|
||||
# zmax = tmax(E, Ec_MeV) * X0_cm
|
||||
# dEdz_max = dEdz(zmax, E) # MeV/cm/nuc
|
||||
# trans_prof = rho_r(r) * dEdz_max # MeV/cm²/nuc
|
||||
# ln, = ax_trans.plot(r, trans_prof, color=color, linewidth=2)
|
||||
# energy_lines.append(ln)
|
||||
|
||||
#aufintegriert über gesamten kristall
|
||||
#physikalisch inkorrekt
|
||||
# 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)) / (2*np.pi*sigma_r**2)
|
||||
# shower_2D = long_prof * trans_prof
|
||||
|
||||
# radial_profile = np.trapz(shower_2D, z, axis=0) # Integration über Tiefe
|
||||
# ln, =ax_trans.plot(r, radial_profile, color=color, linewidth=2)
|
||||
# energy_lines.append(ln)
|
||||
|
||||
#physikalisch korrekt
|
||||
dz = z[1] - z[0]
|
||||
for E, col in zip(energies_MeV, colors):
|
||||
# longitudinales Profil
|
||||
prof_z = dEdz(z, E) # MeV/cm
|
||||
# Integration über z → Energie pro Fläche
|
||||
radial_profile = np.sum(prof_z[:, None] * rho_r(r)[None, :] * dz,axis=0) # MeV/cm²
|
||||
ln, = ax_trans.plot(r, radial_profile, color=col, linewidth=2)
|
||||
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)
|
||||
energy_lines.append(ln)
|
||||
|
||||
# Molière Linien nur positive Linien für Legende
|
||||
# Molière Linien
|
||||
molier_lines = []
|
||||
for r_val, lw in zip(molier_radii, [1.5,2.0]):
|
||||
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)
|
||||
|
||||
ax_trans.set_xlabel("Radius r [cm]", fontsize=fs)
|
||||
|
||||
#ax_trans.set_ylabel(r"$\langle dE/(dz\,dA) \rangle$" + "\n[MeV/cm$^2$/nuc]",fontsize=fs)
|
||||
#ax_trans.set_title("Transverse Profiles at Shower Maximum", fontsize=fs)
|
||||
|
||||
ax_trans.set_ylabel(r"$\int \langle dE/(dz\,dA) \rangle\,dz$" + "\n[MeV/cm$^2$]",fontsize=fs)
|
||||
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)
|
||||
|
|
@ -179,7 +152,7 @@ ax_trans.add_artist(leg2)
|
|||
# Heatmap Legende (Contour)
|
||||
# --------------------------
|
||||
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",
|
||||
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])
|
||||
|
|
|
|||
|
|
@ -7,11 +7,12 @@ general shower theory, different normalisations
|
|||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from scipy.stats import gamma
|
||||
import scienceplots
|
||||
|
||||
# --------------------------
|
||||
# Font size
|
||||
# --------------------------
|
||||
fs = 18
|
||||
fs = 25
|
||||
|
||||
# --------------------------
|
||||
# Material (BGO)
|
||||
|
|
@ -49,14 +50,15 @@ def rho_r(r):
|
|||
z = np.linspace(0, 10, 400)
|
||||
|
||||
plt.figure(figsize=(12, 6))
|
||||
plt.style.use(['science'])
|
||||
|
||||
lines = []
|
||||
for E, lab, col in zip(energies, labels, colors):
|
||||
ln, = plt.plot(z, dEdz(z,E), color=col, linewidth=1.6, label=lab)
|
||||
lines.append(ln)
|
||||
|
||||
plt.xlabel("Depth z [cm]", fontsize=fs)
|
||||
plt.ylabel(r"$\langle dE/dz \rangle$ [MeV/cm/nuc]", fontsize=fs)
|
||||
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 (theory)", fontsize=fs+2)
|
||||
|
||||
plt.xlim(0, z.max())
|
||||
|
|
@ -84,6 +86,7 @@ plt.savefig("plots/BGO_longitudinal_profile_normed.pdf")
|
|||
r = np.linspace(0, 8, 300)
|
||||
|
||||
plt.figure(figsize=(12, 6))
|
||||
plt.style.use(['science'])
|
||||
|
||||
energy_lines = []
|
||||
for E, lab, col in zip(energies, labels, colors):
|
||||
|
|
@ -101,8 +104,8 @@ for E, lab, col in zip(energies, labels, colors):
|
|||
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 [cm]", fontsize=fs)
|
||||
plt.ylabel(r"$\langle dE/(dz\,dA) \rangle$ [MeV/cm$^2$/nuc]", fontsize=fs)
|
||||
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())
|
||||
|
|
@ -143,6 +146,7 @@ plt.savefig("plots/BGO_transverse_profile_normed.pdf")
|
|||
# ============================================================
|
||||
|
||||
plt.figure(figsize=(12, 6))
|
||||
plt.style.use(['science'])
|
||||
dz = z[1]-z[0]
|
||||
energy_lines = []
|
||||
for E, lab, col in zip(energies, labels, colors):
|
||||
|
|
@ -157,8 +161,8 @@ for E, lab, col in zip(energies, labels, colors):
|
|||
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 [cm]", fontsize=fs)
|
||||
plt.ylabel(r"$\int \langle dE/(dz\,dA) \rangle\,dz$" + "\n[MeV/cm$^2$]",fontsize=fs)
|
||||
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)
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -9,6 +9,7 @@ Comparable to analytic shower theory
|
|||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
import scienceplots
|
||||
|
||||
# -------------------------------------------------
|
||||
# Input files
|
||||
|
|
@ -17,28 +18,25 @@ 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",
|
||||
"/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"]
|
||||
|
||||
|
||||
labels = ["100 MeV", "200 MeV", "500 MeV",]
|
||||
colors = ["C0", "C1", "C2"]
|
||||
|
||||
# -------------------------------------------------
|
||||
# Parameters
|
||||
# -------------------------------------------------
|
||||
r_max_cm = 8.0
|
||||
n_bins = 100
|
||||
fs = 18
|
||||
fs = 26
|
||||
|
||||
# -------------------------------------------------
|
||||
# Plot
|
||||
# -------------------------------------------------
|
||||
plt.figure(figsize=(12, 6))
|
||||
plt.style.use(['science'])
|
||||
|
||||
for file, label, color in zip(files, labels, colors):
|
||||
|
||||
|
|
@ -74,13 +72,14 @@ for file, label, color in zip(files, labels, colors):
|
|||
label=label
|
||||
)
|
||||
|
||||
plt.xlabel("Radius r [cm]", fontsize=fs)
|
||||
plt.ylabel(r"$\left\langle \int \frac{dE}{dz\,dA}\,dz \right\rangle$ [MeV/cm$^2$]", fontsize=fs)
|
||||
plt.xlabel("Radius r in cm", fontsize=fs)
|
||||
plt.ylabel(r"$\left\langle \int \frac{\mathrm{d}E}{\mathrm{d}z\,\mathrm{d}A}\,\mathrm{d}z \right\rangle$ in MeV/cm$^2$", fontsize=fs)
|
||||
plt.title("Transverse energy deposition in BGO (z-integrated, Geant4)", fontsize=fs+2)
|
||||
plt.xlim(0, r_max_cm)
|
||||
plt.yscale("log")
|
||||
plt.grid(True, ls="--", lw=0.6)
|
||||
plt.tick_params(labelsize=fs-2)
|
||||
plt.legend(title="Initial Energy")
|
||||
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_integrated.png", dpi=300)
|
||||
|
|
|
|||
|
|
@ -8,6 +8,7 @@ Comparable to analytic shower theory (Plot 2)
|
|||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
import scienceplots
|
||||
|
||||
# -------------------------------------------------
|
||||
# Material & theory parameters (BGO)
|
||||
|
|
@ -43,12 +44,13 @@ n_bins = 80
|
|||
|
||||
dz_slice = 0.5 # cm (slice thickness around shower maximum)
|
||||
|
||||
fs = 18
|
||||
fs = 26
|
||||
|
||||
# -------------------------------------------------
|
||||
# Plot
|
||||
# -------------------------------------------------
|
||||
plt.figure(figsize=(12, 6))
|
||||
plt.style.use(['science'])
|
||||
|
||||
for file, E, label, color in zip(files, energies, labels, colors):
|
||||
|
||||
|
|
@ -100,17 +102,19 @@ for file, E, label, color in zip(files, energies, labels, colors):
|
|||
# -------------------------------------------------
|
||||
# Plot cosmetics
|
||||
# -------------------------------------------------
|
||||
plt.xlabel("Radius r [cm]", fontsize=fs)
|
||||
plt.xlabel("Radius r in cm", fontsize=fs)
|
||||
plt.ylabel(
|
||||
r"$\left\langle \frac{dE}{dz\,dA} \right\rangle$ [MeV/cm$^3$]",
|
||||
r"$\left\langle \frac{\mathrm{d}E}{\mathrm{d}z\,\mathrm{d}A} \right\rangle$ in MeV/cm$^3$",
|
||||
fontsize=fs
|
||||
)
|
||||
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.tick_params(labelsize=fs-2)
|
||||
plt.legend(title="Initial Energy")
|
||||
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)
|
||||
|
|
|
|||
|
|
@ -8,6 +8,7 @@ Step-wise plotting with correct normalization per primary particle
|
|||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
import scienceplots
|
||||
|
||||
# --------------------------
|
||||
# Input files and parameters
|
||||
|
|
@ -23,7 +24,7 @@ files = [
|
|||
energies = [100, 200, 500, 1000, 2000] # MeV
|
||||
colors = ["C0","C1","C2","C3","C4"]
|
||||
|
||||
fs = 14
|
||||
fs = 24
|
||||
|
||||
# --------------------------
|
||||
# Analysis parameters
|
||||
|
|
@ -43,19 +44,20 @@ 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=16, y=0.95)
|
||||
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=12, frameon=True, framealpha=0.85,
|
||||
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 dE/(dz\,dA) \rangle$ [MeV/cm$^2$/prim]",
|
||||
va='center', rotation='vertical', fontsize=16)
|
||||
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
|
||||
|
|
@ -94,15 +96,15 @@ 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.009, None)
|
||||
#ax.set_ylim(0.0009, None)
|
||||
ax.grid(True)
|
||||
ax.tick_params(labelsize=fs)
|
||||
ax.tick_params(labelsize=fs-2)
|
||||
|
||||
if i < n_layers:
|
||||
ax.set_ylabel(f"{z_min}-{z_max} cm", fontsize=fs)
|
||||
else:
|
||||
ax.set_ylabel("Sum", fontsize=fs)
|
||||
ax.set_xlabel("Radius r [cm]", fontsize=fs)
|
||||
ax.set_xlabel("Radius r in cm", fontsize=fs)
|
||||
|
||||
# Legende innerhalb der Achse
|
||||
ax.legend(fontsize=12, frameon=True, framealpha=0.85, facecolor="white")
|
||||
|
|
|
|||
|
|
@ -7,6 +7,7 @@ transversal shower profile theory slicewise
|
|||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from scipy.stats import gamma
|
||||
import scienceplots
|
||||
|
||||
# --------------------------
|
||||
# Material & Parameter
|
||||
|
|
@ -19,6 +20,8 @@ Rm = 2.26
|
|||
energies = [100, 200, 500, 1000, 2000] # MeV
|
||||
colors = ["C0","C1","C2","C3","C4"]
|
||||
|
||||
fs = 24 # globale Schriftgröße
|
||||
|
||||
# --------------------------
|
||||
# Profile
|
||||
# --------------------------
|
||||
|
|
@ -47,19 +50,20 @@ 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=16, y=0.95)
|
||||
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=12, frameon=True, framealpha=0.85,
|
||||
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 dE/(dz\,dA) \rangle$ [MeV/cm$^2$/nuc]", va='center', rotation='vertical', fontsize=16)
|
||||
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
|
||||
|
|
@ -79,13 +83,13 @@ for i, (z_min, z_max) in enumerate(layers + [(0,10)]):
|
|||
ax.set_xlim(0, 8)
|
||||
ax.set_ylim(0, None)
|
||||
ax.grid(True)
|
||||
ax.tick_params(labelsize=14)
|
||||
ax.tick_params(labelsize=fs-2)
|
||||
|
||||
if i < n_layers:
|
||||
ax.set_ylabel(f"{z_min}-{z_max} cm", fontsize=14)
|
||||
ax.set_ylabel(f"{z_min}-{z_max} cm", fontsize=fs)
|
||||
else:
|
||||
ax.set_ylabel("Sum", fontsize=14)
|
||||
ax.set_xlabel("Radius r [cm]", fontsize=14)
|
||||
ax.set_ylabel("Sum", fontsize=fs)
|
||||
ax.set_xlabel("Radius r in cm", fontsize=fs)
|
||||
|
||||
ax.legend(fontsize=12, frameon=True, framealpha=0.85, facecolor="white")
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue