FT_functions.py

@@ -84,69 +84,4 @@ def photon_emission_ava(n,E,m_0,Z,lam1,lam2):
             #print('yes: E=',E[i],";v=",v, "n=", result[i])
         #else: print('no: E=',E[i],";v=",v, "n=")
             
-    return result
+    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)

Frank-Tamm.py

@@ -65,7 +65,6 @@ 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')
@@ -110,65 +109,6 @@ 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)
@@ -254,7 +194,7 @@ fig, ax1 = plt.subplots(figsize=(15, 10))
 plt.style.use(['science'])
 
 # --- linke y-Achse: Cherenkov-Spektrum ---
-ax1.step(
+ax1.plot(
     lam_centers,
     N_lambda,
     lw=2,
@@ -494,7 +434,7 @@ plt.minorticks_on()
 ax2.set_yscale('log')
 #ax1.set_yscale('log')
 
-#plt.savefig("plots/prediction.pdf")
+plt.savefig("plots/prediction.pdf")
 
 
 

n_dependence.py

@@ -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}$ in $\mathrm{cm}^{-1}$',
+    r'$\frac{dN_{\mathrm{photons}}}{dx}\;[\mathrm{cm}^{-1}]$',
     fontsize=fs+2
 )
 
@@ -87,14 +87,6 @@ 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'])
 
@@ -114,18 +106,17 @@ for n_val in n_list:
 
     N_lambda = np.array(N_lambda)
 
-    plt.step(
+    plt.plot(
         lam_centers,
         N_lambda,
         lw=2,
-        color=colors[n_val],
         label=fr'$n={n_val:.2f}$'
     )
 
 # ---------- Plot-Format ----------
-plt.xlabel(r'Wavelength $\lambda$ in nm', fontsize=fs+2)
+plt.xlabel(r'Wavelength $\lambda$ [nm]', fontsize=fs+2)
 plt.ylabel(
-    r'$\frac{dN_{\mathrm{photons}}}{dx}$ in $\mathrm{cm}^{-1}$',
+    r'$\frac{dN_{\mathrm{photons}}}{dx}\;[\mathrm{cm}^{-1}]$',
     fontsize=fs+2
 )
 
@@ -174,6 +165,14 @@ 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 = {}
 
@@ -212,9 +211,9 @@ for n_val in n_list:
     )
 
 # ---------- Linke Achse ----------
-ax1.set_xlabel(r'Wavelength $\lambda$ in nm', fontsize=fs+2)
+ax1.set_xlabel(r'Wavelength $\lambda$ [nm]', fontsize=fs+2)
 ax1.set_ylabel(
-    r'$\frac{\mathrm{d}N_{\mathrm{photons}}}{\mathrm{d}x}$ in $\mathrm{cm}^{-1}$',
+    r'$\frac{\mathrm{d}N_{\mathrm{photons}}}{\mathrm{d}x}\;[\mathrm{cm}^{-1}]$',
     fontsize=fs+2
 )
 ax1.tick_params(axis='both', labelsize=fs)
@@ -236,7 +235,7 @@ ax2.plot(
     label='Transmittance'
 )
 
-ax2.set_ylabel(r'Transmittance in (\%)', fontsize=fs+2)
+ax2.set_ylabel(r'Transmittance (\%)', fontsize=fs+2)
 ax2.tick_params(axis='y', labelsize=fs)
 
 # ---------- Titel ----------
@@ -325,12 +324,7 @@ 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(
+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 ) 
-    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) 
@@ -364,7 +358,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'$E_{\mathrm{thr}}$ in MeV', fontsize=fs+2)
+plt.ylabel(r'Cherenkov threshold energy $E_{\mathrm{thr}}$ [MeV]', fontsize=fs+2)
 
 plt.title(
     r'Cherenkov threshold energy',
@@ -409,8 +403,8 @@ plt.plot(
 
 plt.xlabel(r'Refractive index $n$', fontsize=fs+2)
 plt.ylabel(
-    r'$E_{\mathrm{thr}}$ in GeV',
+    r'Cherenkov threshold energy $E_{\mathrm{thr}}$ [GeV]',
-    fontsize=fs+2
+    fontsize=26
 )
 
 plt.title(

shower_heatmap.py

@@ -9,16 +9,12 @@ 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 in cm")
+    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 in cm")
+    parser.add_argument("--zmax", type=float, default=10.0, help="Max depth [cm]")
     parser.add_argument("--out", default="shower_heatmap.png", help="Output image")
     args = parser.parse_args()
 
@@ -58,16 +54,14 @@ def main():
     )
 
     plt.figure(figsize=(8, 6))
-    im = plt.pcolormesh(R, Z, H_avg, shading="auto", cmap="viridis")
+    im = plt.pcolormesh(R, Z, H_avg, shading="auto", cmap="inferno")
-    cbar = plt.colorbar(im)  # Colorbar erzeugen
+    plt.colorbar(im, label="dE/dz [MeV/cm]")
-    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 in cm", fontsize=fs)
+    plt.xlabel("Radius r [cm]")
-    plt.ylabel("Depth z in cm", fontsize=fs)
+    plt.ylabel("Depth z [cm]")
     plt.ylim(zmax, 0)
-    plt.tick_params(axis='both', which='major', labelsize=fs-2)
+    plt.tick_params(axis='both', which='major', labelsize=fs-5)
-    plt.title("Shower energy deposition 2GeV (Geant4)", fontsize=fs)
+    plt.title("Shower energy deposition (simulation)")
 
     plt.tight_layout()
     plt.savefig(args.out, dpi=300)

shower_long.py

@@ -9,7 +9,6 @@ Comparable to analytic shower theory
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
-import scienceplots 
 
 # -------------------------------------------------
 # Input files
@@ -30,13 +29,12 @@ colors = ["C0", "C1", "C2", "C3", "C4"]
 # -------------------------------------------------
 z_max_cm = 10.0
 n_bins   = 100
-fs = 26
+fs = 18
 
 # -------------------------------------------------
 # Plot
 # -------------------------------------------------
 plt.figure(figsize=(12, 6))
-plt.style.use(['science'])
 
 for file, label, color in zip(files, labels, colors):
 
@@ -65,14 +63,74 @@ for file, label, color in zip(files, labels, colors):
         label=label
     )
 
-plt.xlabel("Depth z in cm", fontsize=fs)
+plt.xlabel("Depth z [cm]", fontsize=fs)
-plt.ylabel(r"$\langle \mathrm{d}E/\mathrm{d}z \rangle$  in MeV/cm", fontsize=fs)
+plt.ylabel(r"$\langle dE/dz \rangle$  [MeV/cm/nuc]", 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", fontsize=fs-4, title_fontsize=fs-2, frameon=True, framealpha=0.85, facecolor="white")
+plt.legend(title="Initial Energy")
 
 plt.tight_layout()
 plt.savefig("plots/G4_longitudinal_profile.png", dpi=300)
-plt.show()
+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()

shower_long_theory.py

@@ -1,13 +1,6 @@
-#!/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
@@ -20,8 +13,6 @@ Rm = 2.26
 energies = [100, 200, 500, 1000, 2000]  # MeV
 colors = ["C0","C1","C2","C3","C4"]
 
-fs = 24  # globale Schriftgröße
-
 # --------------------------
 # Profile
 # --------------------------
@@ -50,13 +41,18 @@ 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=fs+2, y=0.95)
+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))
 
 # Gemeinsames Y-Label
-fig.text(0.02, 0.5, r"$\langle \mathrm{d}E/\mathrm{d}z \rangle$ in MeV/cm", va='center', rotation='vertical', fontsize=fs)
+fig.text(0.02, 0.5, r"$\langle dE/dz \rangle$ [MeV/cm]", va='center', rotation='vertical', fontsize=16)
 
 # --------------------------
 # Plots
@@ -68,9 +64,7 @@ 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")
@@ -78,15 +72,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=fs-2)
+    ax.tick_params(labelsize=14)
     
     if i < n_rings:
-        ax.set_ylabel(f"r={r_min}-{r_max} cm", fontsize=fs)
+        ax.set_ylabel(f"r={r_min}-{r_max} cm", fontsize=14)
     else:
-        ax.set_ylabel("Sum", fontsize=fs)
+        ax.set_ylabel("Sum", fontsize=14)
-        ax.set_xlabel("Depth z in cm", fontsize=fs)
+        ax.set_xlabel("Depth z [cm]", fontsize=14)
     
-    ax.legend(fontsize=fs-10, frameon=True, framealpha=0.85, facecolor="white", loc='upper right')
+    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")

shower_map.py

@@ -1,15 +1,14 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
-Overview: 2D map, transversal and longitudinal shower profile theory
+overwiev: 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
+from scipy.stats import gamma, norm
 from matplotlib.lines import Line2D
-import scienceplots
 
 # --------------------------
 # Funktionen
@@ -43,7 +42,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 = 22  # Schriftgröße
+fs = 16  # Schriftgröße
 
 # Heatmap Konturen in Prozent
 heatmap_levels = [1, 10, 50, 90]
@@ -60,7 +59,6 @@ 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])
@@ -75,13 +73,14 @@ 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], linewidths=2)
+        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 in cm", fontsize=fs)
+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)
@@ -89,7 +88,7 @@ ax_heat.tick_params(labelsize=fs)
 
 # Molière Linien (Heatmap)
 molier_radii = [Rm_cm, 2*Rm_cm]
-molier_widths = [2, 2]
+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)
     ax_heat.axvline(-r_val, color='k', linestyle='-', linewidth=lw)
@@ -103,35 +102,63 @@ 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(r"$\mathrm{d}E/\mathrm{d}z$ in MeV/cm", 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)
+#ax_long.legend(title="Initial Energy", fontsize=fs-2, title_fontsize=fs, loc="upper right", frameon=True)
 
 # --------------------------
-# Transversale Profile (integrated over depth)
+# Transversale Profile (korrekt normiert wie Plot 1.1)
 # --------------------------
 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
-    radial_profile = np.sum(prof_z[:, None] * rho_r(r)[None, :] * dz, axis=0) # MeV/cm²
+    # Integration über z → Energie pro Fläche
-    ln, = ax_trans.step(r, radial_profile, where='mid', color=col, linewidth=2)
+    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)
     energy_lines.append(ln)
 
-# Molière Linien
+# Molière Linien nur positive Linien für Legende
 molier_lines = []
-for r_val, lw in zip(molier_radii, [2,2]):
+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 in cm", fontsize=fs)
+ax_trans.set_xlabel("Radius r [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$",
+#ax_trans.set_ylabel(r"$\langle dE/(dz\,dA) \rangle$" + "\n[MeV/cm$^2$/nuc]",fontsize=fs)
-    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_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)
@@ -152,7 +179,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])

shower_theory.py

@@ -7,12 +7,11 @@ general shower theory, different normalisations
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.stats import gamma
-import scienceplots 
 
 # --------------------------
 # Font size
 # --------------------------
-fs = 25
+fs = 18
 
 # --------------------------
 # Material (BGO)
@@ -50,15 +49,14 @@ 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 in cm", fontsize=fs)
+plt.xlabel("Depth z [cm]", fontsize=fs)
-plt.ylabel(r"$\langle \mathrm{d}E/\mathrm{d}z \rangle$ in MeV/cm", fontsize=fs)
+plt.ylabel(r"$\langle dE/dz \rangle$ [MeV/cm/nuc]", fontsize=fs)
 plt.title("Longitudinal energy deposition in BGO (theory)", fontsize=fs+2)
 
 plt.xlim(0, z.max())
@@ -86,7 +84,6 @@ 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):
@@ -104,8 +101,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 in cm", fontsize=fs)
+plt.xlabel("Radius r [cm]", fontsize=fs)
-plt.ylabel(r"$\langle dE/(\mathrm{d}z\,\mathrm{d}A) \rangle$" + "\nin MeV/cm$^2$", fontsize=fs)
+plt.ylabel(r"$\langle dE/(dz\,dA) \rangle$ [MeV/cm$^2$/nuc]", fontsize=fs)
 plt.title("Transverse energy density at shower maximum (theory)", fontsize=fs+2)
 
 plt.xlim(0, r.max())
@@ -146,7 +143,6 @@ 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):
@@ -161,8 +157,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 in cm", fontsize=fs)
+plt.xlabel("Radius r [cm]", fontsize=fs)
-plt.ylabel(r"$\int \langle dE/(\mathrm{d}z\,\mathrm{d}A) \rangle\,dz$" + "\nin MeV/cm$^2$",fontsize=fs)
+plt.ylabel(r"$\int \langle dE/(dz\,dA) \rangle\,dz$" + "\n[MeV/cm$^2$]",fontsize=fs)
 plt.title("Transverse energy deposition integrated over 10 cm BGO (theory)", fontsize=fs+2)
 
 

shower_trans.py

@@ -9,7 +9,6 @@ Comparable to analytic shower theory
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
-import scienceplots 
 
 # -------------------------------------------------
 # Input files
@@ -18,25 +17,28 @@ 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_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_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 = 26
+fs = 18
 
 # -------------------------------------------------
 # Plot
 # -------------------------------------------------
 plt.figure(figsize=(12, 6))
-plt.style.use(['science'])
 
 for file, label, color in zip(files, labels, colors):
 
@@ -72,14 +74,13 @@ for file, label, color in zip(files, labels, colors):
         label=label
     )
 
-plt.xlabel("Radius r in cm", fontsize=fs)
+plt.xlabel("Radius r [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.ylabel(r"$\left\langle \int \frac{dE}{dz\,dA}\,dz \right\rangle$  [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", fontsize=fs-4, title_fontsize=fs-2, frameon=True, framealpha=0.85, facecolor="white")
+plt.legend(title="Initial Energy")
 
 plt.tight_layout()
 plt.savefig("plots/G4_transverse_profile_integrated.png", dpi=300)

shower_trans_max.py

@@ -8,7 +8,6 @@ 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)
@@ -44,13 +43,12 @@ n_bins   = 80
 
 dz_slice = 0.5   # cm  (slice thickness around shower maximum)
 
-fs = 26
+fs = 18
 
 # -------------------------------------------------
 # Plot
 # -------------------------------------------------
 plt.figure(figsize=(12, 6))
-plt.style.use(['science'])
 
 for file, E, label, color in zip(files, energies, labels, colors):
 
@@ -102,19 +100,17 @@ for file, E, label, color in zip(files, energies, labels, colors):
 # -------------------------------------------------
 # Plot cosmetics
 # -------------------------------------------------
-plt.xlabel("Radius r in cm", fontsize=fs)
+plt.xlabel("Radius r [cm]", fontsize=fs)
 plt.ylabel(
-    r"$\left\langle \frac{\mathrm{d}E}{\mathrm{d}z\,\mathrm{d}A} \right\rangle$  in MeV/cm$^3$",
+    r"$\left\langle \frac{dE}{dz\,dA} \right\rangle$  [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", fontsize=fs-4, title_fontsize=fs-2, frameon=True, framealpha=0.85, facecolor="white")
+plt.legend(title="Initial Energy")
-
 
 plt.tight_layout()
 plt.savefig("plots/G4_transverse_profile_showermax.png", dpi=300)

shower_trans_slices.py

@@ -8,7 +8,6 @@ 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
@@ -24,7 +23,7 @@ files = [
 energies = [100, 200, 500, 1000, 2000]  # MeV
 colors   = ["C0","C1","C2","C3","C4"]
 
-fs = 24
+fs = 14
 
 # --------------------------
 # Analysis parameters
@@ -44,20 +43,19 @@ 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)
+fig.suptitle("Transverse energy deposition in BGO layers (Geant4)", 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]
 axes[0].legend(dummy_lines, [f"{E} MeV" for E in energies],
-               ncol=5, fontsize=fs-10, frameon=True, framealpha=0.85,
+               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 \mathrm{d}E/(\mathrm{d}z\,\mathrm{d}A) \rangle$ in MeV/cm$^2$",
+fig.text(0.02, 0.5, r"$\langle dE/(dz\,dA) \rangle$ [MeV/cm$^2$/prim]",
-         va='center', rotation='vertical', fontsize=fs)
+         va='center', rotation='vertical', fontsize=16)
 
 # --------------------------
 # Layerwise plotting with integrals
@@ -96,15 +94,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.0009, None)
+    ax.set_ylim(0.009, None)
     ax.grid(True)
-    ax.tick_params(labelsize=fs-2)
+    ax.tick_params(labelsize=fs)
 
     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 in cm", fontsize=fs)
+        ax.set_xlabel("Radius r [cm]", fontsize=fs)
 
     # Legende innerhalb der Achse
     ax.legend(fontsize=12, frameon=True, framealpha=0.85, facecolor="white")

shower_trans_theory.py

@@ -7,7 +7,6 @@ transversal shower profile theory slicewise
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.stats import gamma
-import scienceplots 
 
 # --------------------------
 # Material & Parameter
@@ -20,8 +19,6 @@ Rm = 2.26
 energies = [100, 200, 500, 1000, 2000]  # MeV
 colors = ["C0","C1","C2","C3","C4"]
 
-fs = 24  # globale Schriftgröße
-
 # --------------------------
 # Profile
 # --------------------------
@@ -50,20 +47,19 @@ 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)
+fig.suptitle("Transverse energy deposition in BGO layers (theory)", 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=fs-10, frameon=True, framealpha=0.85,
+               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 \mathrm{d}E/(\mathrm{d}z\,\mathrm{d}A) \rangle$ in MeV/cm$^2$/nuc", va='center', rotation='vertical', fontsize=fs)
+fig.text(0.02, 0.5, r"$\langle dE/(dz\,dA) \rangle$ [MeV/cm$^2$/nuc]", va='center', rotation='vertical', fontsize=16)
 
 # --------------------------
 # Plots
@@ -83,13 +79,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=fs-2)
+    ax.tick_params(labelsize=14)
     
     if i < n_layers:
-        ax.set_ylabel(f"{z_min}-{z_max} cm", fontsize=fs)
+        ax.set_ylabel(f"{z_min}-{z_max} cm", fontsize=14)
     else:
-        ax.set_ylabel("Sum", fontsize=fs)
+        ax.set_ylabel("Sum", fontsize=14)
-        ax.set_xlabel("Radius r in cm", fontsize=fs)
+        ax.set_xlabel("Radius r [cm]", fontsize=14)
     
     ax.legend(fontsize=12, frameon=True, framealpha=0.85, facecolor="white")