transverse nad long profiles layerwise correct

shower_heatmap.py

@@ -56,6 +56,7 @@ def main():
     plt.xlabel("Radius r [cm]")
     plt.ylabel("Depth z [cm]")
     plt.ylim(zmax, 0)
+    plt.tick_params(axis='both', which='major', labelsize=fs-5)
     plt.title("Shower energy deposition (simulation)")
 
     plt.tight_layout()

shower_long.py

@@ -54,7 +54,9 @@ for file, label, color in zip(files, labels, colors):
 
 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")
 

shower_long_theory.py

@@ -0,0 +1,88 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.stats import gamma
+
+# --------------------------
+# Material & Parameter
+# --------------------------
+X0 = 1.118
+Ec = 10.5
+b = 0.5
+Rm = 2.26
+
+energies = [100, 200, 500, 1000, 2000]  # MeV
+colors = ["C0","C1","C2","C3","C4"]
+
+# --------------------------
+# Profile
+# --------------------------
+def dEdz(z, E):
+    t = z / X0
+    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)
+
+# --------------------------
+# Gitter
+# --------------------------
+z_grid = np.linspace(0, 10, 400)
+dz = z_grid[1]-z_grid[0]
+r_grid = np.linspace(0, 8, 300)
+dr = r_grid[1]-r_grid[0]
+
+# --------------------------
+# Ringe
+# --------------------------
+rings = [(0,2), (2,4), (4,6), (6,8)]
+n_rings = len(rings)
+
+# --------------------------
+# Figure
+# --------------------------
+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", 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))
+
+# Gemeinsames Y-Label
+fig.text(0.02, 0.5, r"$\langle dE/dz \rangle$ [MeV/cm]", va='center', rotation='vertical', fontsize=16)
+
+# --------------------------
+# Plots
+# --------------------------
+for i, (r_min, r_max) in enumerate(rings + [(0,8)]):
+    ax = axes[i]
+    r_mask = (r_grid >= r_min) & (r_grid < r_max)
+    
+    for E, color in zip(energies, colors):
+        long_prof = dEdz(z_grid, E)
+        radial_profile_2d = long_prof[:, None] * rho_r(r_grid)[None, :]
+        ring_profile = np.trapz(radial_profile_2d[:, r_mask]*2*np.pi*r_grid[r_mask], x=r_grid[r_mask], axis=1)
+        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")
+    
+    ax.set_xlim(0, 10)
+    ax.set_ylim(0, None)
+    ax.grid(True)
+    ax.tick_params(labelsize=14)
+    
+    if i < n_rings:
+        ax.set_ylabel(f"r={r_min}-{r_max} cm", fontsize=14)
+    else:
+        ax.set_ylabel("Sum", fontsize=14)
+        ax.set_xlabel("Depth z [cm]", fontsize=14)
+    
+    ax.legend(fontsize=12, 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")
+plt.savefig("plots/shower_longitudinal_rings_sum.png", dpi=300)
+plt.show()

shower_map.py

@@ -4,7 +4,9 @@ from matplotlib.gridspec import GridSpec
 from scipy.stats import gamma, norm
 from matplotlib.lines import Line2D
 
-# --------- Funktionen ---------
+# --------------------------
+# Funktionen
+# --------------------------
 def tmax(E, Ec):
     return np.log(E / Ec) - 0.5
 
@@ -13,52 +15,53 @@ def longitudinal_profile(t, E, Ec, b=0.5):
     a = b * t_max_val + 1
     return gamma.pdf(t, a, scale=1/b) * E
 
-def transverse_profile(r, Rm, frac=0.9):
-    sigma = Rm / np.sqrt(2 * np.log(1/(1-frac)))
-    return norm.pdf(r, 0, sigma)
+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
 
-# --------- Parameter ---------
+def rho_r(r, Rm=2.26):
+    return np.exp(-r**2/(2*Rm**2)) / (2*np.pi*Rm**2)
+
+# --------------------------
+# Parameter
+# --------------------------
 X0_cm = 1.118
 Ec_MeV = 10.5
 b = 0.5
 Rm_cm = 2.26
 max_depth_cm = 10
+
 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']
 
-# Konturen % der max Energieflussdichte
-heatmap_levels = [1, 10, 50, 90]   
+fs = 16  # Schriftgröße
 
-# --------- Gitter ---------
-# r = np.linspace(-8,8,300)
-# z = np.linspace(0,max_depth_cm,400)
-# Z, R = np.meshgrid(z, r, indexing="ij")
-# --------- Gitter (nur r=0 bis 8) ---------
-r = np.linspace(0, 8, 300)  # statt -8 bis 8
+# Heatmap Konturen in Prozent
+heatmap_levels = [1, 10, 50, 90]
+linestyles = {1:':', 10:'--', 50:'-.', 90:'-'}
+
+# --------------------------
+# Gitter
+# --------------------------
+r = np.linspace(0, 8, 300)
 z = np.linspace(0, max_depth_cm, 400)
 Z, R = np.meshgrid(z, r, indexing="ij")
 
-# --------- Figure und Axes ---------
+# --------------------------
+# Figure und GridSpec
+# --------------------------
 fig = plt.figure(figsize=(14,10))
-fig.suptitle("Electromagnetic shower in BGO (literature)", fontsize=16, y=0.95)
+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)
 
-# --------- Heatmap-Konturen (prozent max) ---------
-linestyles = {1:'-', 10:'--', 50:':', 90:'-.'}
-linestyles = {1:':', 10:'--', 50:'-.', 90:'-'}
-# for E, color in zip(energies_MeV, colors):
-#     t = Z / X0_cm
-#     long_prof = longitudinal_profile(t, E, Ec_MeV, b)
-#     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])
+# --------------------------
+# Heatmap Konturen
+# --------------------------
 for E, color in zip(energies_MeV, colors):
     t = Z / X0_cm
     long_prof = longitudinal_profile(t, E, Ec_MeV, b) / X0_cm
@@ -69,75 +72,85 @@ for E, color in zip(energies_MeV, colors):
     for pct, ls in linestyles.items():
         ax_heat.contour(R, Z, shower_norm, levels=[pct], colors=[color], linestyles=[ls])
 
-
 ax_heat.grid(True, linestyle='--', linewidth=0.5)
-ax_heat.set_ylabel("Depth z [cm]")
-ax_heat.set_ylim(max_depth_cm, 0)  # y=0 oben, z=10 unten
-ax_heat.set_title("2D contourlines")
+ax_heat.set_ylabel("Depth z [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)
+ax_heat.tick_params(labelsize=fs)
 
-
-# Molière Linien
+# Molière Linien (Heatmap)
 molier_radii = [Rm_cm, 2*Rm_cm]
 molier_widths = [1.5, 2.0]
 for r_val, lw in zip(molier_radii, molier_widths):
-    ax_heat.axvline(r_val, color='k', linestyle='-', linewidth=lw, label=f"{r_val} cm")
+    ax_heat.axvline(r_val, color='k', linestyle='-', linewidth=lw)
     ax_heat.axvline(-r_val, color='k', linestyle='-', linewidth=lw)
 
-
-# --------- Longitudinale Profile ---------
-t = np.linspace(0, max_depth_cm/X0_cm, 400) # t in X0
+# --------------------------
+# Longitudinale Profile
+# --------------------------
+t = 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  # MeV/cm
-    ax_long.plot(profile, t*X0_cm, color=color, linewidth=2, label = str(E/1000) + " GeV")
-ax_long.set_title("Longitudinal Profiles")
+    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")
+
+ax_long.set_title("Longitudinal Profiles", fontsize=fs)
+ax_long.set_xlabel("dE/dz [MeV/cm]", fontsize=fs)
 ax_long.grid(True)
-ax_long.set_ylim(max_depth_cm, 0)  # exakt bis 10 cm
-#ax_long.tick_params(labelbottom=False)
-ax_long.set_xlabel("dE/dx [MeV/cm]")
-#ax_long.set_xscale('log')
-ax_long.legend(title="Initial Energy", loc="upper right")
+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 ---------
-# r_trans = np.linspace(-8,8,300)
-# for E, color in zip(energies_MeV, colors):
-#     tmax_X0 = tmax(E, Ec_MeV)
-#     dE_total = longitudinal_profile(tmax_X0, E, Ec_MeV, b)
-#     trans_prof = transverse_profile(np.abs(r_trans), Rm_cm)
-#     trans_prof *= dE_total / np.max(trans_prof)
-#     ax_trans.plot(r_trans, trans_prof, color=color, linewidth=2)
-r_trans = np.linspace(0, 8, 300)  # statt -8 bis 8
+# --------------------------
+# Transversale Profile (korrekt normiert wie Plot 1.1)
+# --------------------------
+r_trans = np.linspace(0, 8, 300)
+
+energy_lines = []
 for E, color in zip(energies_MeV, colors):
-    tmax_X0 = tmax(E, Ec_MeV)
-    dE_total = longitudinal_profile(tmax_X0, E, Ec_MeV, b)/ X0_cm  # MeV/cm
-    trans_prof = transverse_profile(r_trans, Rm_cm)  # kein np.abs() nötig
-    trans_prof *= dE_total / np.max(trans_prof)
-    ax_trans.plot(r_trans, trans_prof, color=color, linewidth=2)
-# molier_radii = [Rm_cm, 2*Rm_cm]
-molier_widths = [1.5, 2.0]
+    zmax = tmax(E, Ec_MeV) * X0_cm
+    dEdz_max = dEdz(zmax, E)  # MeV/cm/nuc
+    trans_prof = rho_r(r_trans) * dEdz_max  # MeV/cm²/nuc
+    ln, = ax_trans.plot(r_trans, trans_prof, color=color, linewidth=2)
+    energy_lines.append(ln)
 
-for r_val, lw in zip(molier_radii, molier_widths):
-    # nur positive Linie bekommt Label für Legende
-    ax_trans.axvline(r_val, color='k', linestyle='-', linewidth=lw, label=f"{r_val:.2f} cm")
+# Molière Linien nur positive Linien für Legende
+molier_lines = []
+for r_val, lw in zip(molier_radii, [1.5,2.0]):
+    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]")
-ax_trans.set_ylabel("dE/dx [MeV/cm]")
-#ax_trans.set_yscale('log')
-ax_trans.set_title("Transverse Profiles")
+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.grid(True)
-ax_trans.legend(title="Molière Radii", loc="upper right")
 ax_trans.set_xlim(0,8)
+ax_trans.tick_params(labelsize=fs)
 
-# --------- Legenden ---------
+# --- 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)
+
+# --------------------------
+# 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", loc='lower right')
+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)
 
-# energy_legend = [Line2D([0],[0], color=col, lw=2) for col in colors]
-# fig.legend(handles=energy_legend, labels=[f"{E} GeV" for E in energies_GeV],
-#            title="Initial Energy", loc='center right', bbox_to_anchor=(0.85,0.2))
-
-plt.tight_layout(rect=[0,0,0.9,0.95])
-plt.savefig("plots/BGO_eshower_math.pdf")
-plt.savefig("plots/BGO_eshower_math.png")
+plt.tight_layout(rect=[0,0,0.95,0.95])
+plt.savefig("plots/shower_map_theory.pdf")
+plt.savefig("plots/shower_map_theory.png", dpi=300)
 plt.show()

shower_map2.py

@@ -0,0 +1,163 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.gridspec import GridSpec
+from scipy.stats import gamma
+from matplotlib.lines import Line2D
+
+# --------------------------
+# 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
+    return gamma.pdf(t, a, scale=1/b) * E  # MeV pro X0
+
+def rho_r(r, Rm=2.26):
+    return np.exp(-r**2/(2*Rm**2)) / (2*np.pi*Rm**2)
+
+# --------------------------
+# Parameter
+# --------------------------
+X0_cm = 1.118
+Ec_MeV = 10.5
+b = 0.5
+Rm_cm = 2.26
+max_depth_cm = 10
+
+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
+
+# Heatmap-Konturen
+heatmap_levels = [1, 10, 50, 90]
+linestyles = {1:':', 10:'--', 50:'-.', 90:'-'}
+
+# --------------------------
+# Gitter
+# --------------------------
+r = np.linspace(0, 8, 300)
+z = np.linspace(0, max_depth_cm, 400)
+Z, R = np.meshgrid(z, r, indexing='ij')
+
+# --------------------------
+# Figure und GridSpec
+# --------------------------
+fig = plt.figure(figsize=(14,10))
+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)
+
+# --------------------------
+# Heatmap + Konturen
+# --------------------------
+energy_lines = []
+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
+
+    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])
+
+    # Dummy-Linie für Initial Energy Legende
+    ln, = ax_heat.plot([], [], color=color, linewidth=2)
+    energy_lines.append(ln)
+
+ax_heat.grid(True, linestyle='--', linewidth=0.5)
+ax_heat.set_ylabel("Depth z [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)
+ax_heat.tick_params(labelsize=fs)
+
+# --------------------------
+# Konturlinien-Legende
+# --------------------------
+line_legend = [Line2D([0],[0], color='k', linestyle=ls, lw=2) for ls in linestyles.values()]
+leg_contour = 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)
+ax_heat.add_artist(leg_contour)
+
+# --------------------------
+# Initial Energy Legende oben rechts (Heatmap)
+# --------------------------
+leg_energy = ax_heat.legend(
+    handles=energy_lines,
+    labels=[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"
+)
+ax_heat.add_artist(leg_energy)
+
+# --------------------------
+# Longitudinale Profile
+# --------------------------
+t = 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)
+ax_long.set_title("Longitudinal Profiles", fontsize=fs)
+ax_long.set_xlabel("dE/dz [MeV/cm]", fontsize=fs)
+ax_long.grid(True)
+ax_long.set_ylim(max_depth_cm, 0)
+ax_long.tick_params(labelsize=fs)
+
+# --------------------------
+# Transversale Profile über gesamte Tiefe
+# --------------------------
+r_trans = np.linspace(0, 8, 300)
+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
+    ax_trans.plot(r_trans, radial_profile, color=color, linewidth=2)
+
+# --------------------------
+# Molier-Radien Linien (Transversal)
+# --------------------------
+molier_radii = [Rm_cm, 2*Rm_cm]
+molier_widths = [1.5, 2.0]
+molier_lines = []
+for r_val, lw in zip(molier_radii, molier_widths):
+    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"$\int \langle dE/(dz\,dA) \rangle$" + "\n[MeV/cm²/nuc]", fontsize=fs)
+ax_trans.set_title("Transverse Profiles integrated over BGO depth", fontsize=fs)
+ax_trans.grid(True)
+ax_trans.set_xlim(0,8)
+ax_trans.tick_params(labelsize=fs)
+
+# --------------------------
+# Legende Molier-Radien (Transversal oben rechts)
+# --------------------------
+leg_molier = 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_trans.add_artist(leg_molier)
+
+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)
+plt.show()

shower_theory.py

@@ -0,0 +1,134 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.stats import gamma
+
+# --------------------------
+# Font size
+# --------------------------
+fs = 16
+
+# --------------------------
+# Material (BGO)
+# --------------------------
+X0 = 1.118      # cm
+Ec = 10.5       # MeV
+b  = 0.5
+Rm = 2.26       # cm
+
+energies = [100, 200, 500, 1000, 2000]  # MeV
+labels   = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
+colors   = ["C0","C1","C2","C3","C4"]
+
+# --------------------------
+# Longitudinal profile
+# --------------------------
+def tmax(E):
+    return np.log(E/Ec) - 0.5
+
+def dEdz(z, E):
+    t = z / X0
+    a = 1 + b * tmax(E)
+    dEdt = E * gamma.pdf(t, a, scale=1/b)
+    return dEdt / X0   # MeV/cm/nuc
+
+# --------------------------
+# Transverse energy density
+# --------------------------
+def rho_r(r):
+    return np.exp(-r**2/(2*Rm**2)) / (2*np.pi*Rm**2)
+
+# ============================================================
+# Plot 1: Longitudinal profile
+# ============================================================
+z = np.linspace(0, 10, 500)
+
+plt.figure(figsize=(12, 6))
+
+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.title("Longitudinal energy deposition in BGO", fontsize=fs)
+
+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.grid(True)
+plt.xticks(fontsize=fs)
+plt.yticks(fontsize=fs)
+plt.tight_layout()
+plt.savefig("plots/BGO_longitudinal_profile_normed.pdf")
+
+# ============================================================
+# Plot 2: Transverse profile at shower maximum
+# ============================================================
+r = np.linspace(0, 8, 400)
+
+plt.figure(figsize=(12, 6))
+
+energy_lines = []
+for E, lab, col in zip(energies, labels, colors):
+    zmax = tmax(E) * X0
+    ln, = plt.plot(
+        r,
+        dEdz(zmax, E) * rho_r(r),
+        color=col,
+        linewidth=1.6,
+        label=lab
+    )
+    energy_lines.append(ln)
+
+# 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 [cm]", fontsize=fs)
+plt.ylabel(r"$\langle dE/(dz\,dA) \rangle$ [MeV/cm$^2$/nuc]", fontsize=fs)
+plt.title("Transverse energy density at shower maximum", fontsize=fs)
+
+plt.xlim(0, r.max())
+plt.ylim(0, None)
+
+# --- Legend 1: Energies ---
+leg1 = plt.legend(
+    handles=energy_lines,
+    title="Initial Energy",
+    fontsize=fs-1,
+    title_fontsize=fs,
+    loc="upper right",
+    frameon=True,
+    framealpha=0.85,
+    facecolor="white"
+)
+
+# --- Legend 2: Geometry ---
+leg2 = plt.legend(
+    handles=[rm_line, rm2_line],
+    fontsize=fs-1,
+    loc="lower right",
+    frameon=True,
+    framealpha=0.85,
+    facecolor="white"
+)
+
+plt.gca().add_artist(leg1)
+
+plt.grid(True)
+plt.xticks(fontsize=fs)
+plt.yticks(fontsize=fs)
+plt.tight_layout()
+plt.savefig("plots/BGO_transverse_profile_normed.pdf")
+plt.show()

shower_trans.py

@@ -45,7 +45,9 @@ for file, label, color in zip(files, labels, colors):
 
 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")
 

shower_trans_slices.py

@@ -0,0 +1,69 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+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()

shower_trans_theory.py

@@ -0,0 +1,89 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.stats import gamma
+
+# --------------------------
+# Material & Parameter
+# --------------------------
+X0 = 1.118
+Ec = 10.5
+b = 0.5
+Rm = 2.26
+
+energies = [100, 200, 500, 1000, 2000]  # MeV
+colors = ["C0","C1","C2","C3","C4"]
+
+# --------------------------
+# Profile
+# --------------------------
+def dEdz(z, E):
+    t = z / X0
+    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)
+
+# --------------------------
+# Gitter
+# --------------------------
+z_grid = np.linspace(0, 10, 400)
+dz = z_grid[1]-z_grid[0]
+r_grid = np.linspace(0, 8, 300)
+dr = r_grid[1]-r_grid[0]
+
+# --------------------------
+# Schichten
+# --------------------------
+layers = [(0,2), (2,4), (4,6), (6,8), (8,10)]
+n_layers = len(layers)
+
+# --------------------------
+# Figure
+# --------------------------
+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", fontsize=16, 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,
+               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)
+
+# --------------------------
+# 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)
+    
+    for E, color in zip(energies, colors):
+        prof = dEdz(z_grid, E)
+        layer_prof = prof[z_mask]
+        radial_profile = np.sum(layer_prof[:, None] * rho_r(r_grid)[None, :] * dz, axis=0)
+        E_layer = np.sum(radial_profile * 2*np.pi*r_grid*dr)
+        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.tick_params(labelsize=14)
+    
+    if i < n_layers:
+        ax.set_ylabel(f"{z_min}-{z_max} cm", fontsize=14)
+    else:
+        ax.set_ylabel("Sum", fontsize=14)
+        ax.set_xlabel("Radius r [cm]", fontsize=14)
+    
+    ax.legend(fontsize=12, frameon=True, framealpha=0.85, facecolor="white")
+
+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)
+plt.show()

sim_map.py

@@ -56,6 +56,8 @@ z_bins = 400
 r_bins = 300
 molier_radii = [2.26, 2*2.26]  # cm
 
+fs=16
+
 # --------- Figure ---------
 fig = plt.figure(figsize=(14,10))
 fig.suptitle("Simulated electromagnetic shower in BGO", fontsize=16, y=0.95)
@@ -172,9 +174,10 @@ for file, label, color in zip(files, energies_labels, colors):
     r_profile_avg[nonzero] /= dr
     ax_trans.step(r_profile_center, r_profile_avg, color=color, lw=2, label=label)
 
-# --------- Heatmap (lokales dE/dz pro Step, OHNE Glättung) ---------
+# Heatmap (lokales dE/dz pro Step, OHNE Glättung) ---------
 
 file=files[2]
+label=energies_labels[2]
 # Grobes, festes Binning: 1 mm = 0.1 cm
 dz = 0.1  # cm
 dr = 0.1  # cm
@@ -217,9 +220,10 @@ im = ax_heat.pcolormesh(
 plt.colorbar(im, ax=ax_heat, label="dE/dz [MeV/cm]")
 ax_heat.set_xlim(0, r_max_cm)
 ax_heat.set_ylim(max_depth_cm, 0)
-ax_heat.set_xlabel("Radius r [cm]")
-ax_heat.set_ylabel("Depth z [cm]")
-ax_heat.set_title(f"Simulated shower {label}")
+ax_heat.set_xlabel("Radius r [cm]",fontsize=fs)
+ax_heat.set_ylabel("Depth z [cm]",fontsize=fs)
+ax_heat.set_title(f"2D shower map {label}",fontsize=fs)
+ax_heat.tick_params(axis='both', which='major', labelsize=fs-3)
 
 
 # --------- Achsen, Titel, Grid ---------
@@ -231,17 +235,19 @@ ax_heat.set_title(f"Simulated shower {label}")
 # ax_heat.grid(True, linestyle='--', linewidth=0.5)
 
 
-ax_long.set_xlabel("dE/dx [MeV/cm]")
+ax_long.set_xlabel("dE/dx [MeV/cm]",fontsize=fs)
 ax_long.set_ylim(max_depth_cm,0)
-ax_long.set_title("Longitudinal Profiles")
+ax_long.set_title("Longitudinal Profiles",fontsize=fs)
 ax_long.grid(True, linestyle='--', linewidth=0.5)
 #ax_long.legend(title="Initial Energy", loc="lower left")
+ax_long.tick_params(axis='both', which='major', labelsize=fs-3)
 
-ax_trans.set_xlabel("Radius r [cm]")
-ax_trans.set_ylabel("dE/dx [MeV/cm]")
-ax_trans.set_title("Transverse Profiles")
+ax_trans.set_xlabel("Radius r [cm]",fontsize=fs)
+ax_trans.set_ylabel("dE/dx [MeV/cm]",fontsize=fs)
+ax_trans.set_title("Transverse Profiles",fontsize=fs)
 ax_trans.grid(True, linestyle='--', linewidth=0.5)
 ax_trans.set_xlim(0,r_max_cm)
+ax_trans.tick_params(axis='both', which='major', labelsize=fs-3)
 ax_trans.legend(title="Initial Energy", loc="upper right")
 
 # --------- Molière Linien ---------