Add theta-angle to Cherenkov plots

FT_functions.py

@@ -84,4 +84,69 @@ 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,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")
 
 
 

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}\;[\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]',
+    r'$E_{\mathrm{thr}}$ in GeV',
-    fontsize=26
+    fontsize=fs+2
 )
 
 plt.title(

shower_heatmap.py

@@ -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("--rmax", type=float, default=8.0, help="Max radius in cm")
-    parser.add_argument("--zmax", type=float, default=10.0, help="Max depth [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")
+    im = plt.pcolormesh(R, Z, H_avg, shading="auto", cmap="viridis")
-    plt.colorbar(im, label="dE/dz [MeV/cm]")
+    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.xlabel("Radius r in cm", fontsize=fs)
-    plt.ylabel("Depth z [cm]")
+    plt.ylabel("Depth z in cm", fontsize=fs)
     plt.ylim(zmax, 0)
-    plt.tick_params(axis='both', which='major', labelsize=fs-5)
+    plt.tick_params(axis='both', which='major', labelsize=fs-2)
-    plt.title("Shower energy deposition (simulation)")
+    plt.title("Shower energy deposition 2GeV (Geant4)", fontsize=fs)
 
     plt.tight_layout()
     plt.savefig(args.out, dpi=300)

shower_long.py

@@ -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.xlabel("Depth z in cm", fontsize=fs)
-plt.ylabel(r"$\langle dE/dz \rangle$  [MeV/cm/nuc]", 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()
+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,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)
+fig.suptitle("Longitudinal energy deposition in BGO by radial ring (theory)", fontsize=fs+2, 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)
+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_ylabel("Sum", fontsize=fs)
-        ax.set_xlabel("Depth z [cm]", fontsize=14)
+        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")

shower_map.py

@@ -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²
-    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)
-    ln, = ax_trans.plot(r, radial_profile, 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_xlabel("Radius r in cm", fontsize=fs)
-
+ax_trans.set_ylabel(
-#ax_trans.set_ylabel(r"$\langle dE/(dz\,dA) \rangle$" + "\n[MeV/cm$^2$/nuc]",fontsize=fs)
+    r"$\int \langle \mathrm{d}E/(\mathrm{d}z\,\mathrm{d}A) \rangle\,\mathrm{d}z$" + "\n in MeV/cm$^2$",
-#ax_trans.set_title("Transverse Profiles at Shower Maximum", fontsize=fs)
+    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)
@@ -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])

shower_theory.py

@@ -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.xlabel("Depth z in cm", fontsize=fs)
-plt.ylabel(r"$\langle dE/dz \rangle$ [MeV/cm/nuc]", 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.xlabel("Radius r in cm", fontsize=fs)
-plt.ylabel(r"$\langle dE/(dz\,dA) \rangle$ [MeV/cm$^2$/nuc]", 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.xlabel("Radius r in cm", fontsize=fs)
-plt.ylabel(r"$\int \langle dE/(dz\,dA) \rangle\,dz$" + "\n[MeV/cm$^2$]",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)
 
 

shower_trans.py

@@ -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_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 = 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.xlabel("Radius r in cm", fontsize=fs)
-plt.ylabel(r"$\left\langle \int \frac{dE}{dz\,dA}\,dz \right\rangle$  [MeV/cm$^2$]", 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)

shower_trans_max.py

@@ -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)

shower_trans_slices.py

@@ -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]",
+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=16)
+         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")

shower_trans_theory.py

@@ -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_ylabel("Sum", fontsize=fs)
-        ax.set_xlabel("Radius r [cm]", fontsize=14)
+        ax.set_xlabel("Radius r in cm", fontsize=fs)
     
     ax.legend(fontsize=12, frameon=True, framealpha=0.85, facecolor="white")