comparinson shower geant4

plot_adriani-combined-e-p-fig.py

@@ -381,7 +381,7 @@ ax.hlines(y=[2.5e-4],xmin=1.0,xmax=11e3,color='tab:blue',alpha=0.2,linewidth=lin
 plt.text(0.91, 0.745, 'HET', fontsize=15, rotation=0, transform=plt.gcf().transFigure, va='center')
 plt.text(0.91, 0.575, 'AMS', fontsize=15, rotation=0, transform=plt.gcf().transFigure, va='center')
 plt.text(0.91, 0.41, 'KET', fontsize=15, rotation=0, transform=plt.gcf().transFigure, va='center')
-plt.text(0.91, 0.24, 'CHAOS', fontsize=15, rotation=0, transform=plt.gcf().transFigure, va='center')
+plt.text(0.91, 0.24, 'ProHEPaM', fontsize=15, rotation=0, transform=plt.gcf().transFigure, va='center')
 
 ax.set_ylim(5e-5,2e0)
 ax.set_xscale('log')
@@ -393,7 +393,7 @@ ax.set_ylabel(r'Differential Flux $\Phi$ in (m$^2$ s sr MeV)$^{-1}$', fontsize=1
 ax.grid(visible=True, which='both', axis='both', ls='--')
 plt.legend(fontsize=12,loc='upper right')
 plt.title('GCR Energy Spectra',size=16)
-plt.savefig('images/adriani-etal-combined-e-p_all_chaos.pdf')
+plt.savefig('plots/adriani-etal-combined-e-p_all_prohepam.pdf')
 
 
 # p_flux = simpson(y1*1000, np.sqrt(x1l*x1h))#, dx = x1h - x1l)

shower.py

@@ -1,221 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from scipy.stats import gamma, norm
-import pandas as pd
-
-# --------------------------
-# Material parameters for BGO
-# --------------------------
-Z_BGO = (4*83 + 3*32 + 12*8)/19  # Approx. effective Z
-X0_cm = 1.118       # Radiation length in cm
-Ec_MeV = 10.5       # Critical energy in MeV
-b = 0.5             # Gamma profile parameter
-dEdx_ion = 9.5  # MeV/cm, approximate ionization loss
-Rm_cm = 2.26         # Approximate Molière radius for BGO
-
-# ----------------------------------------------------------------------------
-# Critical Energy in BGO: Two Definitions
-# --------------------------
-
-# Energy range [MeV]
-E = np.logspace(0, 4, 500)  # 1 MeV to 10 GeV
-
-# --------------------------
-# Losses
-# --------------------------
-dEdx_rad = E / X0_cm
-dEdx_ion_array = np.full_like(E, dEdx_ion)
-dEdx_total = dEdx_rad + dEdx_ion_array
-
-# --------------------------
-# Tsai/Rossi intersection
-# --------------------------
-idx = np.argmin(np.abs(dEdx_rad - dEdx_ion_array))
-Ec_tsai = E[idx]
-dEdx_tsai = dEdx_ion_array[idx]  # Höhe des Schnittpunkts
-
-# --------------------------
-# Rossi definition: dE/dx_ion * (E / dE/dx_rad) = E
-# Equivalent: line parallel to dEdx_rad passing through Ec_rossi
-# --------------------------
-Ec_rossi = dEdx_ion * X0_cm  # MeV
-dEdx_rossi_line = E / Ec_rossi * dEdx_ion  # parallel to radiation line
-
-# --------------------------
-# Plot
-# --------------------------
-# plt.figure(figsize=(8,6))
-
-# # Radiation, ionization, total
-# plt.loglog(E, dEdx_rad, label='Radiation loss dE/dx_rad', linewidth=2, color='blue')
-# plt.loglog(E, dEdx_ion_array, label='Ionization loss dE/dx_ion', linewidth=2, color='orange', linestyle='--')
-# plt.loglog(E, dEdx_total, label='Total loss', linewidth=2, color='k')
-
-# # Tsai/Rossi intersection
-# plt.scatter(Ec_tsai, dEdx_tsai, color='red', s=80, label=f'Tsai/Rossi Ec ≈ {Ec_tsai:.1f} MeV')
-
-# # Rossi definition line (parallel to radiation)
-# plt.loglog(E, dEdx_rossi_line, color='green', linestyle=':', linewidth=2, label=f'Rossi Ec ≈ {Ec_rossi:.1f} MeV')
-
-# plt.xlabel("Energy [MeV]")
-# plt.ylabel("dE/dx [MeV/cm]")
-# plt.title("Critical Energy in BGO: Tsai/Rossi and Rossi Definitions")
-# plt.grid(True, which='both', ls='--', alpha=0.5)
-# plt.legend()
-# plt.tight_layout()
-# plt.show()
-
-
-
-
-# ----------------------------------------------------------------------------
-# Longitudinal energy deposition profile (Gamma function)
-# --------------------------
-
-
-# Primary electron energies [GeV]
-energies_GeV = [0.1, 0.5, 1.0, 2.0]
-energies_MeV = [E*1000 for E in energies_GeV]  # convert to MeV for calculations
-
-# Detector thicknesses [cm]
-thicknesses_cm = [2, 4, 6]
-
-# Font size for all text in plot
-fs = 18
-
-# --------------------------
-# Functions
-# --------------------------
-def tmax(E, Ec):
-    """Calculate shower maximum in radiation lengths"""
-    return np.log(E / Ec) - 0.5
-
-def longitudinal_profile(t, E, Ec, b=0.5):
-    """Longitudinal energy deposition profile (Gamma function)"""
-    t_max_val = tmax(E, Ec)
-    a = b * t_max_val + 1
-    return gamma.pdf(t, a, scale=1/b) * E  # dE/dt in MeV/X0
-
-# --------------------------
-# Calculate table values
-# --------------------------
-rows = []
-
-for E in energies_MeV:
-    t_max_val = tmax(E, Ec_MeV)
-    t_max_cm = t_max_val * X0_cm
-    row = {
-        'tmax [X0]': round(t_max_val,3),
-        'tmax [cm]': round(t_max_cm,3)
-    }
-    for d in thicknesses_cm:
-        thickness_X0 = d / X0_cm
-        row[f'Dicke [X0] ({d}cm)'] = round(thickness_X0,3)
-        row[f'Dicke/tmax ({d}cm)'] = round(thickness_X0 / t_max_val,3)
-    rows.append(row)
-
-# Build a DataFrame with energies as columns (like your LaTeX table)
-table_df = pd.DataFrame(rows, index=[f"{E} GeV" for E in energies_GeV]).T
-# print("\n--- Longitudinal Shower Table ---\n")
-# print(table_df)
-
-# --------------------------
-# Plot longitudinal profiles
-# --------------------------
-t = np.linspace(0, 8, 400)  # Depth in X0
-plt.figure(figsize=(10,6))
-
-for E, E_GeV in zip(energies_MeV, energies_GeV):
-    profile = longitudinal_profile(t, E, Ec_MeV, b)
-    plt.plot(t*X0_cm, profile, label=f"{E_GeV} GeV", linewidth=2)
-
-# Vertical lines for detector thicknesses
-for d in thicknesses_cm:
-    plt.axvline(d, color='k', linestyle='--', alpha=0.5, linewidth=1)
-    plt.text(d+0.05, plt.ylim()[1]*0.9, f"{d} cm", rotation=90, va='top', fontsize=fs)
-
-plt.xlabel("Depth in BGO [cm]", fontsize=fs)
-plt.ylabel("dE/dx [MeV/cm]", fontsize=fs)
-plt.title("Longitudinal Shower Profile in BGO", fontsize=fs)
-plt.legend(fontsize=fs)
-plt.grid(True)
-plt.xticks(fontsize=fs)
-plt.yticks(fontsize=fs)
-plt.tight_layout()
-plt.savefig("plots/BGO_longitudinal_profile.pdf", format='pdf')
-plt.show()
-
-# --------------------------
-# Transverse shower functions (Molière distribution) at shower maximum
-# --------------------------
-def transverse_profile(r, Rm, frac=0.9):
-    """
-    Approximate transverse distribution using Gaussian core of Molière radius
-    frac: fraction of energy within 1 Rm (~90%)
-    """
-    sigma = Rm / np.sqrt(2 * np.log(1/(1-frac)))  # match 90% in Rm
-    return norm.pdf(r, 0, sigma)
-
-r = np.linspace(0, 10, 400)  # transverse distance in cm
-
-
-###PLOT
-plt.figure(figsize=(10,6))
-
-for E, E_GeV in zip(energies_MeV, energies_GeV):
-    tmax_X0 = tmax(E, Ec_MeV)
-    dE_total = longitudinal_profile(tmax_X0, E, Ec_MeV, b)  # peak energy
-    trans_profile = transverse_profile(r, Rm_cm)
-    # normalize to longitudinal peak
-    trans_profile *= dE_total / np.max(trans_profile)
-    plt.plot(r, trans_profile, label=f"{E_GeV} GeV", linewidth=2)
-
-plt.axvline(Rm_cm, color='r', linestyle='--', label="RM="+str(Rm_cm)+"cm")
-plt.axvline(2*Rm_cm, color='r', linestyle=':', label="2RM="+str(2*Rm_cm)+"cm")
-plt.xlabel("Transverse distance [cm]", fontsize=fs)
-plt.ylabel("dE/dx [MeV/cm]", fontsize=fs)
-plt.title("Transverse Shower Profile in BGO at Shower Maximum", fontsize=fs)
-plt.legend(fontsize=fs)
-plt.grid(True)
-plt.xticks(fontsize=fs)
-plt.yticks(fontsize=fs)
-plt.tight_layout()
-plt.show()
-
-#########################################################
-
-# Radius-Achse in cm (symmetrisch um die Schauerachse)
-r = np.linspace(-8, 8, 300)   # 300 radial points
-z = np.linspace(0, 12, 400)   # 400 depth points
-
-# Meshgrid für 2D-Darstellung
-Z, R = np.meshgrid(z, r, indexing="ij")
-
-# Transversales Profil (Gaussian Approximation)
-sigma_r = Rm_cm * (1 + 0.03 * (Z / X0_cm))   # +3% pro X0 (realistische Streuung)
-trans_profile = np.exp(-(R**2) / (2 * sigma_r**2))
-
-t = Z / X0_cm
-E0 = 100     # Beispiel: 1 GeV Elektron
-
-long_profile = longitudinal_profile(t, E0, Ec_MeV, b)
-#long_profile_2D = long_profile[:, None] 
-
-shower_2D = long_profile * trans_profile
-
-plt.figure(figsize=(10, 6))
-plt.imshow(
-    shower_2D,
-    extent=[r.min(), r.max(), z.max(), z.min()],
-    aspect='auto',
-    cmap='inferno',
-    interpolation='bilinear'
-)
-plt.axvline(0, color='black', linestyle='-', linewidth=1.5)
-
-plt.colorbar(label="Energy deposition (arb. units)")
-plt.xlabel("Radius r [cm]")
-plt.ylabel("Depth z [cm]")
-plt.title("2D Heatmap of Electromagnetic Shower in BGO (Rm = 2.26 cm)")
-plt.tight_layout()
-plt.show()

shower_heatmap.py

@@ -1,6 +1,10 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 
+"""
+Heatmap for a chosen file of Geant4 simulation data
+"""
+
 import argparse
 import numpy as np
 import pandas as pd

shower_long.py

@@ -1,11 +1,18 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
+"""
+Longitudinal energy deposition profile from Geant4
+<dE/dz> averaged per primary event
+Comparable to analytic shower theory
+"""
 
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
 
-# --------- Dateien ---------
+# -------------------------------------------------
+# Input files
+# -------------------------------------------------
 files = [
     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_100MeV_0.hits",
     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_200MeV_0.hits",
@@ -17,32 +24,36 @@ files = [
 labels = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
 colors = ["C0", "C1", "C2", "C3", "C4"]
 
-fs=18
-
-# --------- Parameter ---------
+# -------------------------------------------------
+# Parameters
+# -------------------------------------------------
 z_max_cm = 10.0
 n_bins   = 100
+fs = 18
 
-# --------- Plot ---------
+# -------------------------------------------------
+# 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)
+    # number of primary events
+    N_events = df["event"].nunique()
+
+    # z-binning
+    z_edges = np.linspace(0.0, z_max_cm, n_bins + 1)
     dz = z_edges[1] - z_edges[0]
-    z_centers = (z_edges[:-1] + z_edges[1:]) / 2
+    z_centers = 0.5 * (z_edges[:-1] + z_edges[1:])
 
-    # Energie pro z-Bin
+    # energy sum per 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
+    # <dE/dz> [MeV/cm/nuc]
+    profile = E_sum / (N_events * dz)
 
-    # Step-Plot (physikalisch korrekt für Histogramme)
     plt.step(
         z_centers,
         profile,
@@ -53,13 +64,73 @@ 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.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.tick_params(axis='both', which='major', labelsize=fs-3)
-plt.grid(True, ls="--", lw=0.5)
+plt.grid(True, ls="--", lw=0.6)
+plt.tick_params(labelsize=fs-2)
 plt.legend(title="Initial Energy")
 
 plt.tight_layout()
-plt.savefig("plots/shower_longitudinal.png", dpi=300)
+plt.savefig("plots/G4_longitudinal_profile.png", dpi=300)
 plt.show()
+
+
+# # --------- Dateien ---------
+# files = [
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_100MeV_0.hits",
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_200MeV_0.hits",
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_500MeV_0.hits",
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_1GeV_0.hits",
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_2GeV_0.hits",
+# ]
+
+# labels = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
+# colors = ["C0", "C1", "C2", "C3", "C4"]
+
+# fs=18
+
+# # --------- Parameter ---------
+# z_max_cm = 10.0
+# n_bins = 100
+
+# # --------- Plot ---------
+# plt.figure(figsize=(12, 6))
+
+# for file, label, color in zip(files, labels, colors):
+#     df = pd.read_csv(file, sep="\t")
+
+#     z_edges = np.linspace(0, z_max_cm, n_bins + 1)
+#     dz = z_edges[1] - z_edges[0]
+#     z_centers = (z_edges[:-1] + z_edges[1:]) / 2
+
+#     # Energie pro z-Bin
+#     E_sum, _ = np.histogram(df["z"], bins=z_edges, weights=df["edep"])
+#     counts, _ = np.histogram(df["z"], bins=z_edges)
+
+#     # Mittelwert pro Step → dE/dz [MeV/cm]
+#     profile = np.zeros_like(E_sum)
+#     mask = counts > 0
+#     profile[mask] = E_sum[mask] / counts[mask] / dz
+
+#     # Step-Plot (physikalisch korrekt für Histogramme)
+#     plt.step(
+#         z_centers,
+#         profile,
+#         where="mid",
+#         lw=2,
+#         color=color,
+#         label=label
+#     )
+
+# plt.xlabel("Depth z [cm]",fontsize=fs)
+# plt.ylabel("dE/dz [MeV/cm]",fontsize=fs)
+# plt.title("Longitudinal shower profile", fontsize=fs+2)
+# plt.xlim(0, z_max_cm)
+# plt.tick_params(axis='both', which='major', labelsize=fs-3)
+# plt.grid(True, ls="--", lw=0.5)
+# plt.legend(title="Initial Energy")
+
+# plt.tight_layout()
+# plt.savefig("plots/shower_longitudinal.png", dpi=300)
+# plt.show()

shower_long_theory.py

@@ -43,13 +43,13 @@ n_rings = len(rings)
 # --------------------------
 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)
+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))
+# 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)

shower_map.py

@@ -1,3 +1,9 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+overwiev: 2D map, transversal and longitudinal shower profile theory
+"""
+
 import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib.gridspec import GridSpec
@@ -67,6 +73,7 @@ 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():
@@ -104,14 +111,37 @@ ax_long.tick_params(labelsize=fs)
 # --------------------------
 # 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):
-    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)
+
+#lokal im Schauermaximum
+# for E, color in zip(energies_MeV, colors):
+#     zmax = tmax(E, Ec_MeV) * X0_cm
+#     dEdz_max = dEdz(zmax, E)  # MeV/cm/nuc
+#     trans_prof = rho_r(r) * dEdz_max  # MeV/cm²/nuc
+#     ln, = ax_trans.plot(r, trans_prof, color=color, linewidth=2)
+#     energy_lines.append(ln)
+
+#aufintegriert über gesamten kristall
+#physikalisch inkorrekt
+# for E, color in zip(energies_MeV, colors):
+#     t = Z / X0_cm
+#     long_prof = longitudinal_profile(t, E, Ec_MeV, b) / X0_cm
+#     sigma_r = Rm_cm * (1 + 0.03 * t)
+#     trans_prof = np.exp(-(R**2)/(2*sigma_r**2)) / (2*np.pi*sigma_r**2)
+#     shower_2D = long_prof * trans_prof
+
+#     radial_profile = np.trapz(shower_2D, z, axis=0)  # Integration über Tiefe
+#     ln, =ax_trans.plot(r, radial_profile, color=color, linewidth=2)
+#     energy_lines.append(ln)
+
+#physikalisch korrekt
+dz = z[1] - z[0]
+for E, col in zip(energies_MeV, colors):
+    # longitudinales Profil
+    prof_z = dEdz(z, E)                  # MeV/cm
+    # Integration über z → Energie pro Fläche
+    radial_profile = np.sum(prof_z[:, None] * rho_r(r)[None, :] * dz,axis=0) # MeV/cm²
+    ln, = ax_trans.plot(r, radial_profile, color=col, linewidth=2)
     energy_lines.append(ln)
 
 # Molière Linien nur positive Linien für Legende
@@ -122,11 +152,13 @@ for r_val, lw in zip(molier_radii, [1.5,2.0]):
     ax_trans.axvline(-r_val, color='k', linestyle='-', linewidth=lw)
 
 ax_trans.set_xlabel("Radius r [cm]", fontsize=fs)
-ax_trans.set_ylabel(
-    r"$\langle dE/(dz\,dA) \rangle$" + "\n[MeV/cm$^2$/nuc]",
-    fontsize=fs
-)
-ax_trans.set_title("Transverse Profiles at Shower Maximum", fontsize=fs)
+
+#ax_trans.set_ylabel(r"$\langle dE/(dz\,dA) \rangle$" + "\n[MeV/cm$^2$/nuc]",fontsize=fs)
+#ax_trans.set_title("Transverse Profiles at Shower Maximum", fontsize=fs)
+
+ax_trans.set_ylabel(r"$\int \langle dE/(dz\,dA) \rangle\,dz$" + "\n[MeV/cm$^2$]",fontsize=fs)
+ax_trans.set_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)
@@ -151,6 +183,6 @@ ax_heat.legend(line_legend, [f"{pct}%" for pct in linestyles.keys()], title="Con
                fontsize=fs-2, title_fontsize=fs, loc='lower right', frameon=True)
 
 plt.tight_layout(rect=[0,0,0.95,0.95])
-plt.savefig("plots/shower_map_theory.pdf")
-plt.savefig("plots/shower_map_theory.png", dpi=300)
+plt.savefig("plots/shower_map_theory_depth.pdf")
+plt.savefig("plots/shower_map_theory_depth.png", dpi=300)
 plt.show()

shower_map2.py

@@ -1,163 +0,0 @@
-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

@@ -1,3 +1,9 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+general shower theory, different normalisations
+"""
+
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.stats import gamma
@@ -5,7 +11,7 @@ from scipy.stats import gamma
 # --------------------------
 # Font size
 # --------------------------
-fs = 16
+fs = 18
 
 # --------------------------
 # Material (BGO)
@@ -40,7 +46,7 @@ def rho_r(r):
 # ============================================================
 # Plot 1: Longitudinal profile
 # ============================================================
-z = np.linspace(0, 10, 500)
+z = np.linspace(0, 10, 400)
 
 plt.figure(figsize=(12, 6))
 
@@ -51,7 +57,7 @@ for E, lab, col in zip(energies, labels, colors):
 
 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.title("Longitudinal energy deposition in BGO (theory)", fontsize=fs+2)
 
 plt.xlim(0, z.max())
 plt.ylim(0, None)
@@ -75,7 +81,7 @@ plt.savefig("plots/BGO_longitudinal_profile_normed.pdf")
 # ============================================================
 # Plot 2: Transverse profile at shower maximum
 # ============================================================
-r = np.linspace(0, 8, 400)
+r = np.linspace(0, 8, 300)
 
 plt.figure(figsize=(12, 6))
 
@@ -97,7 +103,7 @@ rm2_line = plt.axvline(2*Rm, color='k', linestyle=':',  linewidth=1.6, label=r"$
 
 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.title("Transverse energy density at shower maximum (theory)", fontsize=fs+2)
 
 plt.xlim(0, r.max())
 plt.ylim(0, None)
@@ -131,4 +137,62 @@ plt.xticks(fontsize=fs)
 plt.yticks(fontsize=fs)
 plt.tight_layout()
 plt.savefig("plots/BGO_transverse_profile_normed.pdf")
+
+# ============================================================
+# Plot 3: Transverse profile integrated over entire BGO
+# ============================================================
+
+plt.figure(figsize=(12, 6))
+dz = z[1]-z[0]
+energy_lines = []
+for E, lab, col in zip(energies, labels, colors):
+    # longitudinales Profil
+    prof_z = dEdz(z, E)                  # MeV/cm
+    # Integration über z → Energie pro Fläche
+    radial_profile = np.sum(prof_z[:, None] * rho_r(r)[None, :] * dz,axis=0) # MeV/cm²
+    ln, = plt.plot(r, radial_profile, 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"$\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)
+
+
+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_depth.pdf")
+
 plt.show()

shower_trans.py

@@ -1,56 +1,138 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
+"""
+Transverse energy deposition profile from Geant4
+Integrated over full BGO length (z-integrated)
+Comparable to analytic shower theory
+"""
 
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
 
-# --------- Dateien ---------
+# -------------------------------------------------
+# Input files
+# -------------------------------------------------
 files = [
     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_100MeV_0.hits",
     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_200MeV_0.hits",
     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_500MeV_0.hits",
-    "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_1GeV_0.hits",
-    "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_2GeV_0.hits",
+    #"/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 ---------
+labels = ["100 MeV", "200 MeV", "500 MeV",]
+colors = ["C0", "C1", "C2"]
+
+# -------------------------------------------------
+# Parameters
+# -------------------------------------------------
 r_max_cm = 8.0
 n_bins   = 100
+fs = 18
 
-# --------- Plot ---------
+# -------------------------------------------------
+# Plot
+# -------------------------------------------------
 plt.figure(figsize=(12, 6))
 
 for file, label, color in zip(files, labels, colors):
+
     df = pd.read_csv(file, sep="\t")
+
+    # restrict to BGO radius
     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]
+    # number of primary events
+    N_events = df["event"].nunique()
 
+    # r-binning
+    r_edges = np.linspace(0.0, r_max_cm, n_bins + 1)
+    r_centers = 0.5 * (r_edges[:-1] + r_edges[1:])
+
+    # energy sum per radial bin
     E_sum, _ = np.histogram(df["r"], bins=r_edges, weights=df["edep"])
-    counts, _ = np.histogram(df["r"], bins=r_edges)
 
-    profile = np.zeros_like(E_sum)
-    mask = counts > 0
-    profile[mask] = E_sum[mask] / counts[mask] / dr   # dE/dr [MeV/cm]
+    # ring areas
+    r_in  = r_edges[:-1]
+    r_out = r_edges[1:]
+    A_ring = np.pi * (r_out**2 - r_in**2)
 
-    r_centers = (r_edges[:-1] + r_edges[1:]) / 2
-    plt.step(r_centers, profile, lw=2, color=color, label=label)
+    # <∫ dE / dA dz> [MeV/cm^2]
+    profile = E_sum / (N_events * A_ring)
+
+    plt.step(
+        r_centers,
+        profile,
+        where="mid",
+        lw=2,
+        color=color,
+        label=label
+    )
 
 plt.xlabel("Radius r [cm]", fontsize=fs)
-plt.ylabel("dE/dr [MeV/cm]",fontsize=fs)
-plt.title("Transversal shower profile", fontsize=fs+2)
+plt.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.tick_params(axis='both', which='major', labelsize=fs-5)
-plt.grid(True, ls="--", lw=0.5)
+plt.grid(True, ls="--", lw=0.6)
+plt.tick_params(labelsize=fs-2)
 plt.legend(title="Initial Energy")
 
 plt.tight_layout()
-plt.savefig("plots/shower_transverse.png", dpi=300)
+plt.savefig("plots/G4_transverse_profile_integrated.png", dpi=300)
 plt.show()
+
+
+# # --------- Dateien ---------
+# files = [
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_100MeV_0.hits",
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_200MeV_0.hits",
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_500MeV_0.hits",
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_1GeV_0.hits",
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_2GeV_0.hits",
+# ]
+
+# labels = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
+# colors = ["C0", "C1", "C2", "C3", "C4"]
+
+# fs=18
+
+# # --------- Parameter ---------
+# r_max_cm = 8.0
+# n_bins = 100
+
+# # --------- Plot ---------
+# plt.figure(figsize=(12,6))
+
+# for file, label, color in zip(files, labels, colors):
+#     df = pd.read_csv(file, sep="\t")
+#     df = df[df["r"] <= r_max_cm]
+
+#     r_edges = np.linspace(0, r_max_cm, n_bins + 1)
+#     dr = r_edges[1] - r_edges[0]
+
+#     E_sum, _ = np.histogram(df["r"], bins=r_edges, weights=df["edep"])
+#     counts, _ = np.histogram(df["r"], bins=r_edges)
+
+#     profile = np.zeros_like(E_sum)
+#     mask = counts > 0
+#     profile[mask] = E_sum[mask] / counts[mask] / dr   # dE/dr [MeV/cm]
+
+#     r_centers = (r_edges[:-1] + r_edges[1:]) / 2
+#     plt.step(r_centers, profile, lw=2, color=color, label=label)
+
+# plt.xlabel("Radius r [cm]",fontsize=fs)
+# plt.ylabel("dE/dr [MeV/cm]",fontsize=fs)
+# plt.title("Transversal shower profile", fontsize=fs+2)
+# plt.xlim(0, r_max_cm)
+# plt.tick_params(axis='both', which='major', labelsize=fs-5)
+# plt.grid(True, ls="--", lw=0.5)
+# plt.legend(title="Initial Energy")
+
+# plt.tight_layout()
+# plt.savefig("plots/shower_transverse.png", dpi=300)
+# plt.show()

shower_trans_max.py

@@ -0,0 +1,117 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Transverse energy density at shower maximum from Geant4
+Comparable to analytic shower theory (Plot 2)
+"""
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+# -------------------------------------------------
+# Material & theory parameters (BGO)
+# -------------------------------------------------
+X0 = 1.118   # cm
+Ec = 10.5    # MeV
+b  = 0.5
+
+def tmax(E):
+    """Shower maximum in units of X0"""
+    return np.log(E / Ec) - 0.5
+
+# -------------------------------------------------
+# Input files and energies
+# -------------------------------------------------
+files = [
+    "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_100MeV_0.hits",
+    "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_200MeV_0.hits",
+    "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_500MeV_0.hits",
+    "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_1GeV_0.hits",
+    "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_2GeV_0.hits",
+]
+
+energies = [100, 200, 500, 1000, 2000]  # MeV
+labels   = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
+colors   = ["C0", "C1", "C2", "C3", "C4"]
+
+# -------------------------------------------------
+# Analysis parameters
+# -------------------------------------------------
+r_max_cm = 8.0
+n_bins   = 80
+
+dz_slice = 0.5   # cm  (slice thickness around shower maximum)
+
+fs = 18
+
+# -------------------------------------------------
+# Plot
+# -------------------------------------------------
+plt.figure(figsize=(12, 6))
+
+for file, E, label, color in zip(files, energies, labels, colors):
+
+    df = pd.read_csv(file, sep="\t")
+
+    # Number of primary events
+    N_events = df["event"].nunique()
+
+    # Shower maximum position from theory
+    z_max = tmax(E) * X0
+
+    # Select z-slice around shower maximum
+    df_slice = df[
+        (df["z"] >= z_max - dz_slice/2) &
+        (df["z"] <= z_max + dz_slice/2)
+    ]
+
+    # Restrict to detector radius
+    df_slice = df_slice[df_slice["r"] <= r_max_cm]
+
+    # Radial binning
+    r_edges = np.linspace(0.0, r_max_cm, n_bins + 1)
+    r_centers = 0.5 * (r_edges[:-1] + r_edges[1:])
+
+    # Energy sum per radial bin
+    E_sum, _ = np.histogram(
+        df_slice["r"],
+        bins=r_edges,
+        weights=df_slice["edep"]
+    )
+
+    # Ring areas
+    r_in  = r_edges[:-1]
+    r_out = r_edges[1:]
+    A_ring = np.pi * (r_out**2 - r_in**2)
+
+    # <dE / (dz dA)> [MeV / (cm^3)/nuc]
+    profile = E_sum / (N_events * dz_slice * A_ring)
+
+    plt.step(
+        r_centers,
+        profile,
+        where="mid",
+        lw=2,
+        color=color,
+        label=label
+    )
+
+# -------------------------------------------------
+# Plot cosmetics
+# -------------------------------------------------
+plt.xlabel("Radius r [cm]", fontsize=fs)
+plt.ylabel(
+    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.grid(True, ls="--", lw=0.6)
+plt.tick_params(labelsize=fs-2)
+plt.legend(title="Initial Energy")
+
+plt.tight_layout()
+plt.savefig("plots/G4_transverse_profile_showermax.png", dpi=300)
+plt.show()

shower_trans_slices.py

@@ -1,11 +1,17 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
+"""
+Transverse energy deposition in 2 cm BGO layers from Geant4 simulation
+Step-wise plotting with correct normalization per primary particle
+"""
 
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
 
-# --------- Dateien ---------
+# --------------------------
+# Input files and parameters
+# --------------------------
 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",
@@ -14,56 +20,170 @@ files = [
     "/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"]
+energies = [100, 200, 500, 1000, 2000]  # MeV
 colors   = ["C0","C1","C2","C3","C4"]
 
-fs = 18  # Schriftgröße
+fs = 14
+
+# --------------------------
+# Analysis parameters
+# --------------------------
+z_max_cm = 10.0
 r_max_cm = 8.0
-n_bins = 100
-n_particles = 10000  # Anzahl der simulierten Teilchen
+dz_layer  = 2.0   # cm
+n_bins_r  = 100
 
-# --------- 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
+layers = [(0,2), (2,4), (4,6), (6,8), (8,10)]  # 2 cm slices
+n_layers = len(layers)
 
-# --------- 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)
+r_edges = np.linspace(0, r_max_cm, n_bins_r + 1)
+r_centers = 0.5 * (r_edges[:-1] + r_edges[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):
+# --------------------------
+# 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 (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=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$/prim]",
+         va='center', rotation='vertical', fontsize=16)
+
+# --------------------------
+# Layerwise plotting with integrals
+# --------------------------
+for i, (z_min, z_max) in enumerate(layers + [(0,z_max_cm)]):
+    ax = axes[i]
+
+    for file, E, color in zip(files, energies, colors):
+        df = pd.read_csv(file, sep="\t")
+
+        # Primärereignisse
+        N_events = df["event"].nunique()
+
+        # Slice mask
         df_slice = df[(df["z"] >= z_min) & (df["z"] < z_max)]
+        df_slice = df_slice[df_slice["r"] <= r_max_cm]
+
+        # Histogram radial (Summe der Steps pro Ring)
         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)
+        # Ringflächen
+        r_in  = r_edges[:-1]
+        r_out = r_edges[1:]
+        A_ring = np.pi * (r_out**2 - r_in**2)
 
-# --------- Gemeinsames y-Label ---------
-fig.text(0.02, 0.5, "Energy loss [MeV/cm]", va='center', rotation='vertical', fontsize=fs)
+        # Normierung <dE/(dz dA)> pro Primärteilchen
+        profile = E_sum / (N_events * (z_max - z_min) * A_ring)
 
-# --------- Achsen & Legende ---------
-axes[-1].set_xlabel("Radius r [cm]", fontsize=fs)
-axes[0].legend(title="Initial Energy", loc="upper right", fontsize=fs-4)
+        # Gesamtenergie pro Primärteilchen in dieser Schicht
+        E_layer = E_sum.sum() / N_events
+        linestyle = '--' if i==n_layers else '-'  # Summe gestrichelt
 
-plt.tight_layout(rect=[0.05,0.03,1,0.97])
-plt.savefig("plots/shower_transverse_slices.png", dpi=300)
+        # Step-Plot
+        ax.step(r_centers, profile, where='mid', color=color, linewidth=2,
+                linestyle=linestyle, label=f"{E_layer:.1f} MeV")
+
+    ax.set_xlim(0, r_max_cm)
+    ax.set_yscale('log')
+    ax.set_ylim(0.009, None)
+    ax.grid(True)
+    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 [cm]", fontsize=fs)
+
+    # Legende innerhalb der Achse
+    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/G4_transverse_layers_sum.pdf")
+plt.savefig("plots/G4_transverse_layers_sum.png", dpi=300)
 plt.show()
+
+
+
+
+# #!/usr/bin/env python3
+# # -*- coding: utf-8 -*-
+# """
+# transversal shower profile form Genatz4 simulations slicewise
+# """
+
+# import numpy as np
+# import pandas as pd
+# import matplotlib.pyplot as plt
+
+# # --------- Dateien ---------
+# files = [
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_100MeV_0.hits",
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_200MeV_0.hits",
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_500MeV_0.hits",
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_1GeV_0.hits",
+#     "/home/et189/Geant4/showering/build/out/BGO_shower_e_10-5_2GeV_0.hits",
+# ]
+
+# labels = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
+# colors = ["C0", "C1", "C2", "C3", "C4"]
+
+# fs = 18  # Schriftgröße
+# r_max_cm = 8.0
+# n_bins = 100
+# n_particles = 10000  # Anzahl der simulierten Teilchen
+
+# # --------- Subplots ---------
+# fig, axes = plt.subplots(5, 1, figsize=(8, 12), sharex=True)
+# z_slices = [(0, 2), (2, 4), (4, 6), (6, 8)]
+# dz_slice = 2.0  # Breite der Slice in cm
+
+# # --------- Plotten ---------
+# for file, label, color in zip(files, labels, colors):
+#     df = pd.read_csv(file, sep="\t")
+#     df = df[df["r"] <= r_max_cm]
+
+#     r_edges = np.linspace(0, r_max_cm, n_bins + 1)
+#     dr = r_edges[1] - r_edges[0]
+#     r_centers = (r_edges[:-1] + r_edges[1:]) / 2
+
+#     # --------- Obere 4 Slices ---------
+#     for i, (z_min, z_max) in enumerate(z_slices):
+#         df_slice = df[(df["z"] >= z_min) & (df["z"] < z_max)]
+#         E_sum, _ = np.histogram(df_slice["r"], bins=r_edges, weights=df_slice["edep"])
+#         # Mittelwert pro Teilchen pro cm
+#         profile = E_sum / (n_particles * dz_slice)  
+#         axes[i].step(r_centers, profile, lw=2, color=color, label=label if i==0 else "")
+#         axes[i].set_ylim(0, None)
+#         axes[i].tick_params(axis='both', which='major', labelsize=fs-4)
+#         # kleines y-Label rechts mit Slice
+#         axes[i].text(r_max_cm*1.01, axes[i].get_ylim()[1]*0.9, f"{z_min}-{z_max} cm", rotation=0, va="top", fontsize=fs-6)
+
+# # --------- Unterer Plot: gesamte Summe 0-8cm ---------
+# df_all = df[df["r"] <= r_max_cm]
+# E_sum_total, _ = np.histogram(df_all["r"], bins=r_edges, weights=df_all["edep"])
+# profile_total = E_sum_total / (n_particles * 8.0)  # mittlerer Energieverlust pro cm
+# axes[4].step(r_centers, profile_total, lw=2, color='k')
+# axes[4].set_ylim(0, None)
+# axes[4].tick_params(axis='both', which='major', labelsize=fs-4)
+# axes[4].text(r_max_cm*1.01, axes[4].get_ylim()[1]*0.9, "0-8 cm", rotation=0, va="top", fontsize=fs-6)
+
+# # --------- Gemeinsames y-Label ---------
+# fig.text(0.02, 0.5, "Energy loss [MeV/cm]", va='center', rotation='vertical', fontsize=fs)
+
+# # --------- Achsen & Legende ---------
+# axes[-1].set_xlabel("Radius r [cm]", fontsize=fs)
+# axes[0].legend(title="Initial Energy", loc="upper right", fontsize=fs-4)
+
+# plt.tight_layout(rect=[0.05,0.03,1,0.97])
+# plt.savefig("plots/shower_transverse_slices.png", dpi=300)
+# plt.show()

shower_trans_theory.py

@@ -1,3 +1,9 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+transversal shower profile theory slicewise
+"""
+
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.stats import gamma
@@ -43,7 +49,7 @@ n_layers = len(layers)
 # --------------------------
 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)
+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]

sim_map.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
-Simulationsdaten plotten im Layout wie Literaturplot
+Geant4 Simulationdata map like for literature
 """
 
 import numpy as np