change plotstyle

Frank-Tamm.py

@@ -14,9 +14,18 @@ import scienceplots
 
 from FT_functions import *
 
+import sys
+from pathlib import Path
+
+sys.path.append(str(Path(__file__).resolve().parents[1]))
+import plotstyle
+
+MODE = "paper"
+MODE = "slides"
+plotstyle.set_style(MODE)
+
 ###Frank-Tamm plot 
-plt.figure(figsize=(15, 10))
+plt.figure()
-plt.style.use(['science'])
 
 m_e = 0.511  # in MeV/c^2
 m_p = 938.27  # in MeV/c^2
@@ -82,8 +91,8 @@ plt.plot(E0*J_to_MeV,n_n, color = 'green', lw=2, label = 'helium')
 plt.xscale('log')
 plt.yscale('log')
 
-plt.ylabel(r'$\frac{dN_{\text{photons}}}{dx} \, \text{in cm}^{-1}$', fontsize=fs+2)
+plt.ylabel(r'$\frac{dN_{\text{photons}}}{dx} \, \text{in cm}^{-1}$')
-plt.xlabel(r'$E_{\text{kin}} \, \text{in MeV}$', fontsize=fs+2)
+plt.xlabel(r'$E_{\text{kin}} \, \text{in MeV}$')
 
 plt.xticks(fontsize=fs)
 plt.yticks(fontsize=fs)
@@ -96,24 +105,22 @@ plt.axhline(y=57, color='black', linestyle='--', linewidth=2)
 plt.axhline(y=228, color='green', linestyle='--', linewidth=2)
 
 plt.xlim(0.5,5000000)
+plt.ylim(0.7,400)
 
-plt.title("Number of produced Cherenkov photons in aerogel (300-500nm)",fontsize=fs+3,pad=15)
+plt.title("Number of produced Cherenkov photons in aerogel (300-500nm)")
 
-
+# plt.legend(fontsize=fs-1, loc="lower right",ncol=1,
-#plt.ylim(1,100)
+#     frameon=True, framealpha=0.4,      # leicht transparent
-plt.legend(fontsize=fs-1, loc="lower right",ncol=1,
+#     fancybox=True,
-    frameon=True, framealpha=0.4,      # leicht transparent
+#     #edgecolor='black',   # Rahmenfarbe
-    fancybox=True,
+#     facecolor='white'    # Hintergrundfarbe)
-    #edgecolor='black',   # Rahmenfarbe
+# )
-    facecolor='white'    # Hintergrundfarbe)
+plt.legend(loc="lower right")
-)
+plotstyle.savefig("frank-tamm", category="Cherenkov")
-plt.savefig("plots/frank-tamm")
-#plt.show()
 
 
 ##################### ANGLE #########################
-plt.figure(figsize=(15, 10))
+plt.figure()
-plt.style.use(['science'])
 
 # Cherenkov angles
 theta_e  = cherenkov_theta_deg(1.05, E0, m_e)
@@ -149,22 +156,18 @@ plt.text(1, theta_max_105, rf'n=1.05: ${theta_max_105:.1f}^\circ$', color='black
 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.xlabel(r'$E_{\mathrm{kin}}\;\text{in MeV}$')
-plt.ylabel(r'Cherenkov angle $\theta_C$ in deg', fontsize=fs+2)
+plt.ylabel(r'Cherenkov angle $\theta_C$ in deg')
 
-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.title('Cherenkov angle in aerogel')
 
-plt.legend(fontsize=fs-2, loc='lower right',
+plt.legend(loc='lower right')
-           frameon=True, framealpha=0.4, fancybox=True)
+plotstyle.savefig("cherenkov_angle_deg", category="Cherenkov")
-
-plt.savefig("plots/cherenkov_angle_deg")
 
 
 ###############################################################
@@ -200,27 +203,18 @@ for l1, l2 in zip(Nlam[:-1], Nlam[1:]):
 N_lambda = np.array(N_lambda)
 
 # ---------- Plot ----------
-plt.figure(figsize=(15, 10))
+plt.figure()
-plt.style.use(['science'])
 
 plt.step(lam_centers, N_lambda, color = 'blue', where='mid',lw=2)
 #plt.scatter(lam_centers, N_lambda, s=40, color = 'blue')
 
-plt.xticks(fontsize=fs)
-plt.yticks(fontsize=fs)
-plt.grid(alpha=0.3)
-
 plt.title(
-    r'Cherenkov-Spektrum in aerogel for an $100\,\mathrm{MeV}$ electron',
+    r'Cherenkov-Spektrum in aerogel for an $100\,\mathrm{MeV}$ electron')
-    fontsize=fs+3
+plt.ylabel(r'$\frac{\mathrm{d}N_{\text{photons}}}{\mathrm{d}x}\,\mathrm{cm}^{-1}$')
-)
+plt.xlabel(r'Wavelength $\lambda$ in nm')
-plt.ylabel(r'$\frac{\mathrm{d}N_{\text{photons}}}{\mathrm{d}x}\,\mathrm{cm}^{-1}$',
-           fontsize=fs + 2)
-plt.xlabel(r'Wavelength $\lambda$ in nm', fontsize=fs+2)
 
 plt.xlim(200, 800)
-plt.tight_layout()
+plotstyle.savefig("spectrum", category="Cherenkov")
-plt.savefig("plots/spectrum.pdf")
 
 
 ##################  TRANSMITTENZ  ###############
@@ -250,8 +244,7 @@ T = T_percent / 100.0
 T_percent_centers = 0.5 * (T_percent[:-1] + T_percent[1:])
 T_centers = T_percent_centers / 100.0
 
-fig, ax1 = plt.subplots(figsize=(15, 10))
+fig, ax1 = plt.subplots()
-plt.style.use(['science'])
 
 # --- linke y-Achse: Cherenkov-Spektrum ---
 ax1.step(
@@ -269,10 +262,9 @@ ax1.set_ylabel(
     fontsize=fs+2, color='blue'
 )
 ax1.tick_params(axis='y', labelcolor='blue')
-ax1.tick_params(axis='both', labelsize=fs)
+#ax1.tick_params(axis='both', labelsize=fs)
-ax1.grid(alpha=0.3)
 ax1.tick_params(axis='y', labelcolor='blue', labelsize=fs)
-ax1.grid(alpha=0.3)
+
 
 # --- rechte y-Achse: Transmittanz ---
 ax2 = ax1.twinx()
@@ -294,14 +286,12 @@ ax1.plot(
     label=r'$\left(\mathrm{d}N/\mathrm{d}x\right)\times T$'
 )
 
-ax2.set_ylabel(r'Transmittance (\%)', fontsize=fs+2, color='red')
+ax2.set_ylabel(r'Transmittance (\%)', color='red')
-ax2.tick_params(axis='y', labelcolor='red', labelsize=fs)
+ax2.tick_params(axis='y', labelcolor='red')
 
 # --- Titel & Limits ---
 ax1.set_title(
-    r'Effective Cherenkov photon yield in aerogel for a $100\,\mathrm{MeV}$ electron',
+    r'Effective Cherenkov photon yield in aerogel for a $100\,\mathrm{MeV}$ electron')
-    fontsize=fs+5
-)
 
 ax1.set_xlim(260, 700)
 ax1.set_ylim(0, 6.5)
@@ -325,24 +315,12 @@ labels  = [l_dndx[0], l_T[0], l_dndx[1]]
 ax1.legend(
     handles,
     labels,
-    fontsize=fs,
+    loc='center right'
-    loc='center right',
-    frameon=True,   
-    fancybox=True,
-    #edgecolor='black',   # Rahmenfarbe
-    facecolor='white'    # Hintergrundfarbe
 )
 
-plt.xticks(fontsize=fs)
-plt.yticks(fontsize=fs)
-
-plt.grid(True, which='major', alpha=0.9, linewidth=1.0)
-plt.grid(True, which='minor', alpha=0.8, linewidth=0.5)
-
 #plt.minorticks_on()
 
-plt.savefig("plots/Transspec.pdf")
+plotstyle.savefig("Transspec", category="Cherenkov")
-
 
 #################      QUANTEN EFFIZIENZ        #####################
 QE_percent = np.array([
@@ -361,33 +339,10 @@ QE_percent = np.array([
 
 mask = QE_percent[:-1] >= 0
 
-plt.figure(figsize=(8, 16))
-plt.style.use(['science'])
-
-#plt.step(Nlam[:-1][mask],QE_percent[:-1][mask],where='mid',lw=2)
-plt.plot(lam_centers[mask], QE_percent[:-1][mask], 'x')
-
-plt.xlabel("Wavelength $\lambda$ in nm", fontsize=fs+2)
-plt.ylabel("Quantum Efficiency in \%", fontsize=fs+2)
-plt.title("Photocathode Quantum Efficiency", fontsize=fs+5)
-plt.xlim(200, 800)
-plt.yscale('log')
-
-plt.xticks(fontsize=fs)
-plt.yticks(fontsize=fs)
-
-plt.grid(True, which='major', alpha=0.9, linewidth=1.0)
-plt.grid(True, which='minor', alpha=0.8, linewidth=0.5)
-
-plt.minorticks_on()
-plt.savefig("plots/QE.pdf")
-
-
 ####### PREDICTION ##########
 
 
-fig, ax1 = plt.subplots(figsize=(15, 10))
+fig, ax1 = plt.subplots()
-plt.style.use(['science'])
 
 # --- linke y-Achse: Cherenkov-Spektrum ---
 ax1.plot(
@@ -399,16 +354,13 @@ ax1.plot(
 )
 
 
-ax1.set_xlabel(r'Wavelength $\lambda$ in nm', fontsize=fs+2)
+ax1.set_xlabel(r'Wavelength $\lambda$ in nm')
 ax1.set_ylabel(
-    r'$\frac{\mathrm{d}N_{\text{photons}}}{\mathrm{d}x}\,\mathrm{cm}^{-1}$',
+    r'$\frac{\mathrm{d}N_{\text{photons}}}{\mathrm{d}x}\,\mathrm{cm}^{-1}$', color='blue'
-    fontsize=fs+2, color='blue'
 )
 ax1.tick_params(axis='y', labelcolor='blue')
-ax1.tick_params(axis='both', labelsize=fs)
+ax1.tick_params(axis='both')
-ax1.grid(alpha=0.3)
+ax1.tick_params(axis='y', labelcolor='blue')
-ax1.tick_params(axis='y', labelcolor='blue', labelsize=fs)
-ax1.grid(alpha=0.3)
 
 # --- rechte y-Achse: Transmittanz ---
 ax2 = ax1.twinx()
@@ -453,14 +405,12 @@ ax1.fill_between(
 )
 
 
-ax2.set_ylabel(r'Transmittance (\%)', fontsize=fs+2, color='red')
+ax2.set_ylabel(r'Transmittance (\%)', color='red')
-ax2.tick_params(axis='y', labelcolor='red', labelsize=fs)
+ax2.tick_params(axis='y', labelcolor='red')
 
 # --- Titel & Limits ---
 ax1.set_title(
-    r'Effective Cherenkov photon yield in aerogel for a $100\,\mathrm{MeV}$ electron',
+    r'Effective Cherenkov photon yield in aerogel for a $100\,\mathrm{MeV}$ electron')
-    fontsize=fs+5
-)
 
 ax1.set_xlim(260, 700)
 #ax1.set_ylim(0, 6.5)
@@ -475,15 +425,30 @@ labels  = [l_dndx[0], l_T[0], l_T[1], l_dndx[1], l_dndx[2], l_dndx[3]]
 ax1.legend(
     handles,
     labels,
-    fontsize=fs,
     loc='upper right',
-    frameon=True,
-    fancybox=True,
-    #framealpha=0.9,      # leicht transparent
-    edgecolor='black',   # Rahmenfarbe
-    facecolor='white'    # Hintergrundfarbe
 )
 
+plt.minorticks_on()
+ax2.set_yscale('log')
+#ax1.set_yscale('log')
+
+plotstyle.savefig("prediction_wrong", category="Cherenkov")
+
+
+####### QE ################
+
+plt.figure(figsize=(8, 16))
+
+#plt.step(Nlam[:-1][mask],QE_percent[:-1][mask],where='mid',lw=2)
+plt.plot(lam_centers[mask], QE_percent[:-1][mask], 'x')
+
+plt.xlabel("Wavelength $\lambda$ in nm", fontsize=fs+2)
+plt.ylabel("Quantum Efficiency in \%", fontsize=fs+2)
+plt.title("Photocathode Quantum Efficiency", fontsize=fs+5)
+plt.xlim(200, 800)
+plt.yscale('log')
+plt.ylim(0.008,50)
+
 plt.xticks(fontsize=fs)
 plt.yticks(fontsize=fs)
 
@@ -491,10 +456,8 @@ plt.grid(True, which='major', alpha=0.9, linewidth=1.0)
 plt.grid(True, which='minor', alpha=0.8, linewidth=0.5)
 
 plt.minorticks_on()
-ax2.set_yscale('log')
+#plt.savefig("plots/QE.pdf")
-#ax1.set_yscale('log')
+plotstyle.savefig("QE", category="Cherenkov")
-
-#plt.savefig("plots/prediction.pdf")
 
 
 

n_dependence.py

@@ -6,7 +6,15 @@ import scienceplots
 
 from FT_functions import *
 
-fs=40
+import sys
+from pathlib import Path
+
+sys.path.append(str(Path(__file__).resolve().parents[1]))
+import plotstyle
+
+MODE = "paper"
+#MODE = "slides"
+plotstyle.set_style(MODE)
 
 #### Dependence on n #########
 
@@ -34,6 +42,10 @@ z = z_e
 lam1 = 300
 lam2 = 500
 
+# =================================================
+# Photonenzahl vs n
+# =================================================
+
 # Brechungsindex-Bereich
 n_values = np.linspace(1.04, 1.10, 100)
 
@@ -52,34 +64,18 @@ for n_val in n_values:
 
 N_total = np.array(N_total)
 
-# ---------- Plot ----------
+plt.figure()
-plt.figure(figsize=(15, 10))
+plt.plot(n_values, N_total, lw=2, color='blue', label='electron')
-plt.style.use(['science'])
+plt.xlabel(r'Refractive index $n$')
-
+plt.ylabel(r'$\frac{dN_{\mathrm{photons}}}{dx}$ in $\mathrm{cm}^{-1}$')
-plt.plot(n_values, N_total, lw=2, color='blue')
+plt.title(r'Cherenkov photon yield ($100\,\mathrm{MeV}$ electron)')
-
+plt.legend()
-plt.xlabel(r'Refractive index $n$', fontsize=fs+2)
+plotstyle.savefig("photons_n", category="Cherenkov")
-plt.ylabel(
-    r'$\frac{dN_{\mathrm{photons}}}{dx}$ in $\mathrm{cm}^{-1}$',
-    fontsize=fs+2
-)
-
-plt.title(
-    r'Cherenkov photon yield'
-    '\n$100\,\mathrm{MeV}$ electron, 300–500 nm',
-    fontsize=fs+3
-)
-
-plt.grid(alpha=0.3)
-plt.xticks(fontsize=fs)
-plt.yticks(fontsize=fs)
-
-plt.tight_layout()
-plt.savefig("plots/photons_n.pdf")
 
 
-#########   PHOTON EMISSSION PRO WELLENLÄNGE   ###############################
+# =================================================
-
+# Spektrum vs n
+# =================================================
 
 # ---------- Wellenlängenbins ----------
 Nlam = np.arange(200, 810, 10)   # nm
@@ -95,8 +91,7 @@ colors = {
     1.09: 'tab:red'
 }
 
-plt.figure(figsize=(15, 10))
+plt.figure()
-plt.style.use(['science'])
 
 for n_val in n_list:
     N_lambda = []
@@ -122,29 +117,17 @@ for n_val in n_list:
         label=fr'$n={n_val:.2f}$'
     )
 
-# ---------- Plot-Format ----------
+plt.xlabel(r'Wavelength $\lambda$ in nm')
-plt.xlabel(r'Wavelength $\lambda$ in nm', fontsize=fs+2)
+plt.ylabel(r'$\frac{dN_{\mathrm{photons}}}{dx}$ in $\mathrm{cm}^{-1}$')
-plt.ylabel(
+plt.title(r'Cherenkov spectrum in aerogel ($100\,\mathrm{MeV}$ electron)')
-    r'$\frac{dN_{\mathrm{photons}}}{dx}$ in $\mathrm{cm}^{-1}$',
-    fontsize=fs+2
-)
-
-plt.title(
-    'Cherenkov spectrum in aerogel ($100\,\mathrm{MeV}$ electron)',
-    fontsize=fs+3
-)
-
 plt.xlim(200, 800)
-plt.grid(alpha=0.3)
+plt.legend()
-plt.legend(fontsize=fs-2, frameon=True)
+plotstyle.savefig("spectrum_n", category="Cherenkov")
 
-plt.xticks(fontsize=fs)
+# =================================================
-plt.yticks(fontsize=fs)
+# Spektrum × Transmittanz
+# =================================================
 
-plt.tight_layout()
-plt.savefig("plots/spectrum_n.pdf")
-
-################# SPEKTRUM UND TRANSMITTENZ ###########
 T_percent = np.array([
     0, 0, 0, 0.2, 1.1, 2.5, 4.5, # 200–260
     7.1, 10.4, 14.4, 18.8,    # 270–300
@@ -171,9 +154,7 @@ T = T_percent / 100.0
 T_percent_centers = 0.5 * (T_percent[:-1] + T_percent[1:])
 T_centers = T_percent_centers / 100.0
 
-fig, ax1 = plt.subplots(figsize=(15, 10))
+fig, ax1 = plt.subplots()
-plt.style.use(['science'])
-
 # ---------- Cherenkov-Spektren ----------
 spectra = {}
 
@@ -211,62 +192,23 @@ for n_val in n_list:
         #label=fr'$n={n_val:.2f},\ (\mathrm{{d}}N/\mathrm{{d}}x)\times T$'
     )
 
-# ---------- Linke Achse ----------
+ax1.set_xlabel(r'Wavelength $\lambda$ in nm')
-ax1.set_xlabel(r'Wavelength $\lambda$ in nm', fontsize=fs+2)
+ax1.set_ylabel(r'$\frac{dN_{\mathrm{photons}}}{dx}$ in $\mathrm{cm}^{-1}$')
-ax1.set_ylabel(
-    r'$\frac{\mathrm{d}N_{\mathrm{photons}}}{\mathrm{d}x}$ in $\mathrm{cm}^{-1}$',
-    fontsize=fs+2
-)
-ax1.tick_params(axis='both', labelsize=fs)
-ax1.grid(True, which='major', alpha=0.8, linewidth=1.0)
-ax1.grid(True, which='minor', alpha=0.5, linewidth=0.5)
-ax1.minorticks_on()
-
 ax1.set_xlim(260, 700)
 ax1.set_ylim(0, 11)
-# ---------- Rechte Achse: Transmittanz ----------
+
 ax2 = ax1.twinx()
+ax2.plot(lam_centers[mask], T_percent[:-1][mask], lw=2, linestyle='--', color='black', label='Transmittance')
+ax2.set_ylabel(r'Transmittance in (\%)')
 
-ax2.plot(
-    lam_centers[mask],
-    T_percent[:-1][mask],
-    linestyle='--',
-    color='black',
-    lw=2,
-    label='Transmittance'
-)
-
-ax2.set_ylabel(r'Transmittance in (\%)', fontsize=fs+2)
-ax2.tick_params(axis='y', labelsize=fs)
-
-# ---------- Titel ----------
-ax1.set_title(
-    r'Cherenkov spectrum and effective photon yield in aerogel'
-    '\n$100\,\mathrm{MeV}$ electron',
-    fontsize=fs+4
-)
-
-# ---------- Gemeinsame Legende ----------
 h1, l1 = ax1.get_legend_handles_labels()
 h2, l2 = ax2.get_legend_handles_labels()
+ax1.legend(h1 + h2, l1 + l2, loc='center right')
+plotstyle.savefig("Transspec_multi_n", category="Cherenkov")
 
-ax1.legend(
+# =================================================
-    h1 + h2,
+# Effektiver Photonenzahl
-    l1 + l2,
+# =================================================
-    fontsize=fs-2,
-    loc='upper right',
-    frameon=True,
-    fancybox=True,
-    framealpha=0.85,
-    facecolor='white'
-)
-
-plt.tight_layout()
-plt.savefig("plots/Transspec_multi_n.pdf")
-
-
-
-################### Vergleichbare Phototnen mit Transmittenz ####
 
 lam1 = 300
 lam2 = 500
@@ -321,30 +263,16 @@ for n_val in n_values:
 
 N_eff_total = np.array(N_eff_total)
 
-plt.figure(figsize=(15, 10)) 
+plt.figure()
-plt.style.use(['science']) 
+plt.plot(n_values, N_eff_total, lw=2, color='black')
-plt.plot( n_values, N_eff_total, lw=2, color='black' ) 
+plt.xlabel(r'Refractive index $n$')
-plt.xlabel(r'Refractive index $n$', fontsize=fs+2) 
+plt.ylabel(r'Effective Cherenkov photon yield ')
-plt.ylabel(
+plt.title(r'Cherenkov photon yield (300–500 nm) vs. refractive index')
-    r'Effective Cherenkov photon yield'
+plotstyle.savefig("photons_n_effective", category="Cherenkov")
-    '\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) 
-plt.minorticks_on() 
-plt.xticks(fontsize=fs) 
-plt.yticks(fontsize=fs) 
-plt.tight_layout() 
-plt.savefig("plots/photons_n_effective.pdf")
 
-
+# =================================================
-
+# Schwellenenergie für alle Teilchen
-
+# =================================================
-
-##### Schwellenenergien ########
 
 E_thr_e  = threshold_energy(n_values, m_e)  * J_to_MeV
 E_thr_p  = threshold_energy(n_values, m_p)  * J_to_MeV
@@ -352,40 +280,21 @@ E_thr_mu = threshold_energy(n_values, m_mu) * J_to_MeV
 E_thr_a  = threshold_energy(n_values, m_a)  * J_to_MeV
 
 
-# ---------- Plot ----------
+plt.figure()
-plt.figure(figsize=(15, 10))
+plt.plot(n_values, E_thr_e, color='blue', lw=2, label='electron')
-plt.style.use(['science'])
+plt.plot(n_values, E_thr_p, color='red', lw=2, label='proton')
-
+plt.plot(n_values, E_thr_a, color='green', lw=2, label='helium')
-plt.plot(n_values, E_thr_e,  lw=2, color='blue',  label='electron')
+plt.plot(n_values, E_thr_mu, color=(0.5, 0.3, 0.8), lw=2, label='muon')
-plt.plot(n_values, E_thr_p,  lw=2, color='red',   label='proton')
-plt.plot(n_values, E_thr_a,  lw=2, color='green', label='helium')
-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$')
+plt.ylabel(r'$E_{\mathrm{thr}}$ in MeV')
+plt.title(r'Cherenkov threshold energy')
+plt.legend()
+plotstyle.savefig("threshold_energy_n", category="Cherenkov")
 
-plt.xlabel(r'Refractive index $n$', fontsize=fs+2)
+# =================================================
-plt.ylabel(r'$E_{\mathrm{thr}}$ in MeV', fontsize=fs+2)
+# Schwellenenergie Proton (Zoom)
-
+# =================================================
-plt.title(
-    r'Cherenkov threshold energy',
-    fontsize=30,
-    pad=15
-)
-
-plt.grid(alpha=0.3)
-plt.legend(
-    fontsize=fs-2,
-    frameon=True,
-    framealpha=0.4,
-    fancybox=True
-)
-
-plt.xticks(fontsize=fs)
-plt.yticks(fontsize=fs)
-plt.minorticks_on()
-plt.tight_layout()
-plt.savefig("plots/threshold_energy_n.pdf")
-
 
 # ---------- n-Bereich ----------
 n_values_p = np.linspace(1.05, 1.10, 200)
@@ -393,47 +302,13 @@ n_values_p = np.linspace(1.05, 1.10, 200)
 # ---------- Schwellenenergie Proton ----------
 E_thr_p = threshold_energy(n_values_p, m_p) * J_to_MeV
 
-# ---------- Plot ----------
+plt.figure()
-plt.figure(figsize=(15, 10))
+plt.plot(n_values_p, E_thr_p/1000, color='red', lw=2, label='proton')
-plt.style.use(['science'])
+plt.xlabel(r'Refractive index $n$')
-
+plt.ylabel(r'$E_{\mathrm{thr}}$ in GeV')
-plt.plot(
+plt.title(r'Cherenkov threshold energy for protons')
-    n_values_p,
+plt.legend()
-    E_thr_p/1000,
+plotstyle.savefig("threshold_energy_proton", category="Cherenkov")
-    lw=2,
-    color='red',
-    label='proton'
-)
-
-#plt.yscale('log')
-
-plt.xlabel(r'Refractive index $n$', fontsize=fs+2)
-plt.ylabel(
-    r'$E_{\mathrm{thr}}$ in GeV',
-    fontsize=fs+2
-)
-
-plt.title(
-    r'Cherenkov threshold energy for protons'
-    '\n$1.05 \le n \le 1.10$',
-    fontsize=fs+5,
-    pad=15
-)
-
-plt.grid(alpha=0.3)
-plt.legend(
-    fontsize=fs-2,
-    frameon=True,
-    framealpha=0.4,
-    fancybox=True
-)
-
-plt.xticks(fontsize=fs)
-plt.yticks(fontsize=fs)
-plt.grid(which='minor', linestyle='--', linewidth=0.5, alpha=0.5)
-plt.minorticks_on()
-plt.tight_layout()
-plt.savefig("plots/threshold_energy_proton.pdf")
 
 
 plt.show()

plot_adriani-combined-e-p-fig.py

@@ -2,7 +2,16 @@ import matplotlib.pyplot as plt
 import numpy as np
 from scipy.integrate import simpson, quad
 import scienceplots
-plt.style.use(['science'])
+import sys
+from pathlib import Path
+
+sys.path.append(str(Path(__file__).resolve().parents[1]))
+import plotstyle
+
+MODE = "paper"
+#MODE = "slides"
+plotstyle.set_style(MODE)
+
 color_cycle = plt.rcParams['axes.prop_cycle'].by_key()['color']
 
 def FFS(T=np.logspace(-2,2,100), a=0.4): # Force Field Solution
@@ -28,7 +37,7 @@ x2l, x2h, y2 = np.loadtxt('datafiles/adriani-etal-2015-electron-spectra.dat', sk
 
 #print(np.sqrt(x1l*x1h), x1l, x1h)
 
-figsize = (8,5)
+#figsize = (8,5)
 linewidth = 12
 # fig, ax = plt.subplots(figsize=(6,3.5))
 
@@ -352,7 +361,7 @@ linewidth = 12
 
 
 
-fig, ax = plt.subplots(figsize=figsize)
+fig, ax = plt.subplots()
 
 #energy ranges (xmin and xmax are not in data coordinates!!!)
 xlow, xhigh = 1e-2, 1e2
@@ -360,6 +369,8 @@ def scale(x):
     return (x-xlow)/(xhigh-xlow)
 ax.set_xlim(xlow,xhigh)
 
+linewidth=22
+
 ax.plot(np.sqrt(x1l*x1h), y1, c='k', ls='none', label='Protons p$^+$', marker = 'o',zorder=3)
 ax.plot(np.sqrt(x2l*x2h), y2/1000, c='k', ls='none', label='Electrons e$^-$', marker = 's',zorder=3)
 ax.hlines(y=[4e-1],xmin=10e-3,xmax=107e-3,color='tab:orange',alpha=0.5,linewidth=linewidth,label='Range p$^+$')
@@ -378,22 +389,25 @@ ax.hlines(y=[2.5e-4],xmin=1.0,xmax=11e3,color='tab:blue',alpha=0.2,linewidth=lin
 # ax.plot(T_a/1000, sp*1.e4, c='b', ls='-', label='FFS')
 # ax.plot(T_a/1000, se*1.e4, c='b', ls='-')
 
-plt.text(0.91, 0.745, 'HET', fontsize=15, rotation=0, transform=plt.gcf().transFigure, va='center')
+plt.text(0.91, 0.745, 'HET', fontsize=20, 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.575, 'AMS', fontsize=20, 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.41, 'KET', fontsize=20, 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')
+plt.text(0.91, 0.24, 'ProHEPaM', fontsize=20, rotation=0, transform=plt.gcf().transFigure, va='center')
 
 ax.set_ylim(5e-5,2e0)
 ax.set_xscale('log')
 ax.set_yscale('log')
-ax.set_xlabel(r'Kinetic Energy $E_{\rm kin}$ in GeV',fontsize=14)
+ax.set_xlabel(r'Kinetic Energy $E_{\rm kin}$ in GeV')
-ax.set_ylabel(r'Differential Flux $\Phi$ in (m$^2$ s sr MeV)$^{-1}$', fontsize=14)
+ax.set_ylabel(r'Differential Flux $\Phi$ in (m$^2$ s sr MeV)$^{-1}$')
 # ax.set_xlabel(r'$E_{\rm kin}$ in GeV',fontsize=14)
 # ax.set_ylabel(r'd$J$/d$E$ in (m$^2$ s sr MeV)$^{-1}$', fontsize=14)
+
 ax.grid(visible=True, which='both', axis='both', ls='--')
-plt.legend(fontsize=12,loc='upper right')
+plt.legend(loc='upper right')
-plt.title('GCR Energy Spectra',size=16)
+plt.title('GCR Energy Spectra')
-plt.savefig('plots/adriani-etal-combined-e-p_all_prohepam.pdf')
+plotstyle.savefig('adriani-etal-combined-e-p_all_prohepam', category="other")
+
+plt.show()
 
 
 # p_flux = simpson(y1*1000, np.sqrt(x1l*x1h))#, dx = x1h - x1l)

shower_long.py

@@ -69,10 +69,11 @@ plt.xlabel("Depth z in cm", fontsize=fs)
 plt.ylabel(r"$\langle \mathrm{d}E/\mathrm{d}z \rangle$  in MeV/cm", fontsize=fs)
 plt.title("Longitudinal energy deposition in BGO (Geant4)", fontsize=fs+2)
 plt.xlim(0, z_max_cm)
-plt.grid(True, ls="--", lw=0.6)
+plt.grid(True)
 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.tight_layout()
 plt.savefig("plots/G4_longitudinal_profile.png", dpi=300)
-plt.show()
+plt.savefig("plots/G4_longitudinal_profile.pdf")
+plt.show()

shower_long_comparison.py

@@ -0,0 +1,202 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Longitudinal shower profile comparison
+Theory (analytic, dashed) vs Geant4 simulation (steps)
+With residuals panel
+"""
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from scipy.stats import gamma
+from matplotlib.lines import Line2D
+
+from pathlib import Path
+import sys
+sys.path.append(str(Path(__file__).resolve().parents[1]))
+import plotstyle
+
+import scienceplots  # optional, falls noch andere Style-Features gebraucht werden
+
+# =====================================================
+# Plot Style
+# =====================================================
+MODE = "slides"  # "paper" oder "slides"
+plotstyle.set_style(MODE)
+
+# -------------------------------------------------
+# Input files (Geant4)
+# -------------------------------------------------
+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"]
+
+# -------------------------------------------------
+# Material (BGO)
+# -------------------------------------------------
+X0 = 1.118  # cm
+Ec = 10.5   # MeV
+b  = 0.5
+
+# -------------------------------------------------
+# Plot parameters
+# -------------------------------------------------
+z_max_cm = 10.0
+n_bins   = 100
+
+# -------------------------------------------------
+# Theory functions
+# -------------------------------------------------
+def tmax(E):
+    return np.log(E / Ec) - 0.5
+
+def dEdz_theory(z, E):
+    t = z / X0
+    a = 1 + b * tmax(E)
+    return E * gamma.pdf(t, a, scale=1 / b) / X0  # MeV/cm/nuc
+
+# =====================================================
+# Figure layout
+# =====================================================
+fig, (ax, ax_res) = plt.subplots(
+    2, 1,
+    figsize=(plt.rcParams["figure.figsize"][0], plt.rcParams["figure.figsize"][1]),
+    sharex=True,
+    gridspec_kw={"height_ratios": [3, 1], "hspace": 0.08}
+)
+
+energy_handles = []
+
+# -------------------------------------------------
+# Loop energies
+# -------------------------------------------------
+for file, E, label, color in zip(files, energies, labels, colors):
+
+    # -------------------------
+    # Geant4 simulation
+    # -------------------------
+    df = pd.read_csv(file, sep="\t")
+    N_events = df["event"].nunique()
+
+    z_edges = np.linspace(0.0, z_max_cm, n_bins + 1)
+    dz = z_edges[1] - z_edges[0]
+    z_centers = 0.5 * (z_edges[:-1] + z_edges[1:])
+
+    E_sum, _ = np.histogram(df["z"], bins=z_edges, weights=df["edep"])
+    profile_sim = E_sum / (N_events * dz)
+
+    ax.step(
+        z_centers,
+        profile_sim,
+        where="mid",
+        lw=plt.rcParams["lines.linewidth"],
+        color=color,
+        alpha=0.9
+    )
+
+    # -------------------------
+    # Theory
+    # -------------------------
+    z_th = np.linspace(0, z_max_cm, 600)
+    profile_th = dEdz_theory(z_th, E)
+
+    ln, = ax.plot(
+        z_th,
+        profile_th,
+        color=color,
+        linewidth=plt.rcParams["lines.linewidth"],
+        linestyle="--"
+    )
+
+    energy_handles.append(Line2D([0], [0], color=color, lw=plt.rcParams["lines.linewidth"], linestyle="-"))
+
+    # -------------------------
+    # Residuals
+    # -------------------------
+    th_interp = np.interp(z_centers, z_th, profile_th)
+    mask = th_interp > 1e-6
+
+    residuals = np.zeros_like(profile_sim)
+    residuals[mask] = profile_sim[mask]-th_interp[mask]
+
+    ax_res.step(
+        z_centers,
+        residuals,
+        color=color,
+        linewidth=plt.rcParams["lines.linewidth"]/1.5
+    )
+
+    # schwarze Linie bei y=0
+    ax_res.axhline(0.0, color="black", linestyle=":", linewidth=plt.rcParams["lines.linewidth"]/1.5)
+
+# -------------------------------------------------
+# Main axis
+# -------------------------------------------------
+ax.set_ylabel(r"$\langle \mathrm{d}E/\mathrm{d}z \rangle$ in MeV/cm", fontsize=plt.rcParams["axes.labelsize"])
+ax.set_title("Longitudinal Energy Deposition in BGO", fontsize=plt.rcParams["axes.titlesize"], pad=plt.rcParams.get("axes.titlepad", 15))
+ax.set_xlim(0, z_max_cm)
+ax.grid(True)
+ax.tick_params(labelsize=plt.rcParams["xtick.labelsize"])
+
+# -------------------------------------------------
+# Residual axis
+# -------------------------------------------------
+ax_res.set_xlabel("Depth z in cm", fontsize=plt.rcParams["axes.labelsize"])
+ax_res.set_ylabel(r"$\mathrm{G4sim} - \mathrm{Theory}$", fontsize=plt.rcParams["ytick.labelsize"])
+ax_res.axhline(0.0, color="black", linestyle=":", linewidth=plt.rcParams["lines.linewidth"]/1.5)
+ax_res.grid(True, which="both", linestyle="--", linewidth=0.4)
+ax_res.tick_params(labelsize=plt.rcParams["xtick.labelsize"])
+ax_res.set_xlim(0, z_max_cm)
+ax_res.set_ylim(-20,20)
+
+# -------------------------------------------------
+# Legends
+# -------------------------------------------------
+# Energy legend (colors only)
+leg1 = ax.legend(
+    handles=energy_handles,
+    labels=labels,
+    title="Initial Energy",
+    fontsize=plt.rcParams["legend.fontsize"]-2,
+    title_fontsize=plt.rcParams["legend.fontsize"],
+    loc="upper right",
+    frameon=True,
+    framealpha=plt.rcParams.get("legend.framealpha", 0.85),
+    facecolor="white"
+)
+
+# Style legend
+style_handles = [
+    Line2D([0], [0], color="black", linewidth=plt.rcParams["lines.linewidth"], linestyle="--"),
+    Line2D([0], [0], color="black", linewidth=plt.rcParams["lines.linewidth"], linestyle="-"),
+]
+
+style_labels = ["Theory", "Geant4 simulation"]
+
+leg2 = ax.legend(
+    handles=style_handles,
+    labels=style_labels,
+    fontsize=plt.rcParams["legend.fontsize"]-3,
+    loc="upper left",
+    frameon=True,
+    framealpha=plt.rcParams.get("legend.framealpha", 0.85),
+    facecolor="white"
+)
+
+ax.add_artist(leg1)
+
+# -------------------------------------------------
+# Save
+# -------------------------------------------------
+plt.tight_layout()
+plotstyle.savefig("shower_long_comparison", category="BGO")
+plt.show()

shower_map.py

@@ -1,7 +1,8 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
-Overview: 2D map, transversal and longitudinal shower profile theory
+2D map, transversal and longitudinal shower profile theory (BGO, NKG model)
+Paper-style version, fully cleaned
 """
 
 import numpy as np
@@ -9,153 +10,150 @@ import matplotlib.pyplot as plt
 from matplotlib.gridspec import GridSpec
 from scipy.stats import gamma
 from matplotlib.lines import Line2D
-import scienceplots
+from pathlib import Path
+import sys
+import scienceplots 
+
+# Plotstyle importieren
+sys.path.append(str(Path(__file__).resolve().parents[1]))
+import plotstyle
+
+# ============================================================
+# Paper-style setzen
+# ============================================================
+plotstyle.set_style("paper")
 
 # --------------------------
-# Funktionen
+# Physikalische 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 * tmax(E, Ec) + 1
-    a = b * t_max_val + 1
     return gamma.pdf(t, a, scale=1/b) * E
 
 def dEdz(z, E, X0=1.118, b=0.5, Ec=10.5):
-    t = z / X0
     a = 1 + b * (np.log(E/Ec) - 0.5)
-    return E * gamma.pdf(t, a, scale=1/b) / X0   # MeV/cm/nuc
+    return E * gamma.pdf(z/X0, a, scale=1/b) / X0
 
-def rho_r(r, Rm=2.26):
+def rho_r(r, Rm_cm=2.26, s=1.0):
-    return np.exp(-r**2/(2*Rm**2)) / (2*np.pi*Rm**2)
+    x = r / Rm_cm
+    core = np.zeros_like(r)
+    mask = r > 0
+    core[mask] = (x[mask]**(s-2)) * ((1 + x[mask])**(s-4.5))
+    if r[0] == 0:
+        core[0] = ((1e-9 / Rm_cm)**(s-2)) * ((1 + 1e-9 / Rm_cm)**(s-4.5))
+    norm = np.trapz(2*np.pi*r*core, r)
+    return core / norm
 
 # --------------------------
-# Parameter
+# Parameter (BGO)
 # --------------------------
 X0_cm = 1.118
 Ec_MeV = 10.5
-b = 0.5
 Rm_cm = 2.26
 max_depth_cm = 10
+b = 0.5
 
 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
+heatmap_levels = [0.1, 1, 10]
-
+linestyles = [":", "--", "-"]
-# Heatmap Konturen in Prozent
-heatmap_levels = [1, 10, 50, 90]
-linestyles = {1:':', 10:'--', 50:'-.', 90:'-'}
 
 # --------------------------
 # Gitter
 # --------------------------
-r = np.linspace(0, 8, 300)
+r_max = 2.5
+r = np.linspace(0, r_max, 300)
 z = np.linspace(0, max_depth_cm, 400)
 Z, R = np.meshgrid(z, r, indexing="ij")
 
 # --------------------------
-# Figure und GridSpec
+# Figure und Layout
 # --------------------------
-fig = plt.figure(figsize=(14,10))
+fig = plt.figure()
-plt.style.use(['science'])
+gs = GridSpec(2, 2, width_ratios=[4, 1], height_ratios=[4, 1],
-fig.suptitle("Electromagnetic shower in BGO (theory)", fontsize=fs+2, y=0.97)
+              hspace=0.3, wspace=0.1)
-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_heat = fig.add_subplot(gs[0, 0])
-ax_long = fig.add_subplot(gs[0,1], sharey=ax_heat)
+ax_long = fig.add_subplot(gs[0, 1], sharey=ax_heat)
-ax_trans = fig.add_subplot(gs[1,0], sharex=ax_heat)
+ax_trans = fig.add_subplot(gs[1, 0], sharex=ax_heat)
 
 # --------------------------
-# Heatmap Konturen
+# 2D-Konturen
 # --------------------------
 for E, color in zip(energies_MeV, colors):
-    t = Z / X0_cm
+    long_prof = dEdz(z, E)
-    long_prof = longitudinal_profile(t, E, Ec_MeV, b) / X0_cm
+    shower_2D = np.zeros((len(z), len(r)))
-    sigma_r = Rm_cm * (1 + 0.03 * t)
+    rho_rz = rho_r(r)
-    trans_prof = np.exp(-(R**2)/(2*sigma_r**2))
+    rho_rz[0] = rho_rz[1]
-    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.grid(True, linestyle='--', linewidth=0.5)
+    for i in range(len(z)):
-ax_heat.set_ylabel("Depth z in cm", fontsize=fs)
+        shower_2D[i, :] = long_prof[i] * rho_rz
+
+    levels = np.array(heatmap_levels)/100 * np.max(shower_2D)
+    ax_heat.contour(R, Z, shower_2D,
+                    levels=levels,
+                    colors=[color]*len(levels),
+                    linestyles=linestyles)
+
+# Achsen
+ax_heat.set_ylabel("Depth z in cm")
 ax_heat.set_ylim(max_depth_cm, 0)
-ax_heat.set_xlim(0,8)
+ax_heat.set_xlim(0, r_max)
-ax_heat.set_title("2D contourlines", fontsize=fs)
+ax_heat.set_title("2D Energy density contours (normalized)")
-ax_heat.tick_params(labelsize=fs)
 
-# Molière Linien (Heatmap)
+# Molière-Linie
-molier_radii = [Rm_cm, 2*Rm_cm]
+molier_radii = [Rm_cm]
-molier_widths = [2, 2]
+for r_val in molier_radii:
-for r_val, lw in zip(molier_radii, molier_widths):
+    ax_heat.axvline(r_val, color='k', linestyle='-')
-    ax_heat.axvline(r_val, color='k', linestyle='-', linewidth=lw)
-    ax_heat.axvline(-r_val, color='k', linestyle='-', linewidth=lw)
 
 # --------------------------
-# Longitudinale Profile
+# Longitudinalprofile
 # --------------------------
-t = np.linspace(0, max_depth_cm/X0_cm, 400)
+t_plot = 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
+    profile = longitudinal_profile(t_plot, E, Ec_MeV, b) / X0_cm
-    ax_long.plot(profile, t*X0_cm, color=color, linewidth=2, label=f"{E/1000:.1f} GeV")
+    ax_long.plot(profile, t_plot*X0_cm, color=color, label=f"{E/1000:.1f} GeV")
 
-ax_long.set_title("Longitudinal Profiles", fontsize=fs)
+ax_long.set_title("Longitudinal Profiles")
-ax_long.set_xlabel(r"$\mathrm{d}E/\mathrm{d}z$ in MeV/cm", fontsize=fs)
+ax_long.set_xlabel(r"$\mathrm{d}E/\mathrm{d}z$ in MeV/cm")
-ax_long.grid(True)
 ax_long.set_ylim(max_depth_cm, 0)
-ax_long.tick_params(labelsize=fs)
+ax_long.set_xlim(left=0)
 
 # --------------------------
-# Transversale Profile (integrated over depth)
+# Transversales Profil
 # --------------------------
+n_bins = 80
+r_edges = np.linspace(0, r_max, n_bins+1)
+r_centers = 0.5*(r_edges[:-1]+r_edges[1:])
+z_grid = np.linspace(0, max_depth_cm, 400)
+dz = z_grid[1]-z_grid[0]
+
 energy_lines = []
-dz = z[1] - z[0]
 for E, col in zip(energies_MeV, colors):
-    prof_z = dEdz(z, E)                  # MeV/cm
+    prof_z = dEdz(z_grid, E)
-    radial_profile = np.sum(prof_z[:, None] * rho_r(r)[None, :] * dz, axis=0) # MeV/cm²
+    rho = rho_r(r_centers)
-    ln, = ax_trans.step(r, radial_profile, where='mid', color=col, linewidth=2)
+    profile = np.sum(prof_z[:, None] * rho[None, :] * dz, axis=0)
+    ln, = ax_trans.plot(r_centers, profile, color=col, label=f"{E/1000:.1f} GeV")
     energy_lines.append(ln)
 
-# Molière Linien
+for r_val in molier_radii:
-molier_lines = []
+    ax_trans.axvline(r_val, color='k', linestyle='-')
-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 in cm", fontsize=fs)
+ax_trans.set_yscale("log")
-ax_trans.set_ylabel(
+ax_trans.set_ylim(0.5, None)
-    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 integrated over 10cm")
-    fontsize=fs
+ax_trans.set_xlabel("Radius r in cm")
-)
+ax_trans.set_ylabel(r"$\int \langle \mathrm{d}E/(\mathrm{d}z\,\mathrm{d}A) \rangle\, \mathrm{d}z$")
-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)
-
-# --- Legenden ---
-leg1 = ax_heat.legend(energy_lines, [f"{E/1000:.1f} GeV" for E in energies_MeV],
-                       title="Initial Energy", fontsize=fs-2, title_fontsize=fs,
-                       loc="upper right", frameon=True, framealpha=0.85, facecolor="white")
-
-leg2 = ax_trans.legend(molier_lines, [f"{r_val:.2f} cm" for r_val in molier_radii],
-                       title="Molière Radii", fontsize=fs-2, title_fontsize=fs,
-                       loc="upper right", frameon=True, framealpha=0.85, facecolor="white")
-
-ax_heat.add_artist(leg1)
-ax_trans.add_artist(leg2)
 
 # --------------------------
-# Heatmap Legende (Contour)
+# Legenden
 # --------------------------
-line_legend = [Line2D([0],[0], color='k', linestyle=ls, lw=2) for ls in linestyles.values()]
+ax_heat.legend(energy_lines, [f"{E/1000:.1f} GeV" for E in energies_MeV], title="Initial Energy", loc='upper right')
-ax_heat.legend(line_legend, [f"{pct}\%" for pct in linestyles.keys()], title="Contour % of max",
+plt.tight_layout()
-               fontsize=fs-2, title_fontsize=fs, loc='lower right', frameon=True)
+plotstyle.savefig("shower_map_theory_depth", category="BGO")
-
-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

@@ -39,10 +39,18 @@ def dEdz(z, E):
     return dEdt / X0   # MeV/cm/nuc
 
 # --------------------------
-# Transverse energy density
+# NKG-like transverse profile; safe at r=0
 # --------------------------
+s = 1.0
 def rho_r(r):
-    return np.exp(-r**2/(2*Rm**2)) / (2*np.pi*Rm**2)
+    x = r / Rm
+    core = np.zeros_like(r)
+    mask = r > 0
+    core[mask] = (x[mask]**(s-2)) * ((1 + x[mask])**(s-4.5))
+    if r[0] == 0:
+        core[0] = ((1e-9 / Rm)**(s-2)) * ((1 + 1e-9 / Rm)**(s-4.5))
+    norm = np.trapz(2*np.pi*r*core, r)
+    return core / norm
 
 # ============================================================
 # Plot 1: Longitudinal profile
@@ -64,15 +72,8 @@ plt.title("Longitudinal energy deposition in BGO (theory)", fontsize=fs+2)
 plt.xlim(0, z.max())
 plt.ylim(0, None)
 
-leg = plt.legend(
+plt.legend(handles=lines, title="Initial Energy", fontsize=fs-1, title_fontsize=fs,
-    handles=lines,
+           frameon=True, framealpha=0.85, facecolor="white")
-    title="Initial Energy",
-    fontsize=fs-1,
-    title_fontsize=fs,
-    frameon=True,
-    framealpha=0.85,
-    facecolor="white"
-)
 
 plt.grid(True)
 plt.xticks(fontsize=fs)
@@ -83,7 +84,7 @@ plt.savefig("plots/BGO_longitudinal_profile_normed.pdf")
 # ============================================================
 # Plot 2: Transverse profile at shower maximum
 # ============================================================
-r = np.linspace(0, 8, 300)
+r = np.linspace(0, 2.5, 300)  # nur bis 2.5 cm
 
 plt.figure(figsize=(12, 6))
 plt.style.use(['science'])
@@ -102,14 +103,14 @@ for E, lab, col in zip(energies, labels, colors):
 
 # 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 in cm", fontsize=fs)
 plt.ylabel(r"$\langle dE/(\mathrm{d}z\,\mathrm{d}A) \rangle$" + "\nin MeV/cm$^2$", fontsize=fs)
 plt.title("Transverse energy density at shower maximum (theory)", fontsize=fs+2)
 
-plt.xlim(0, r.max())
+plt.xlim(0, 2.5)
-plt.ylim(0, None)
+plt.yscale("log")  # logarithmische y-Achse
+plt.ylim(0.1, 100000)  # minimale y für log
 
 # --- Legend 1: Energies ---
 leg1 = plt.legend(
@@ -125,21 +126,20 @@ leg1 = plt.legend(
 
 # --- Legend 2: Geometry ---
 leg2 = plt.legend(
-    handles=[rm_line, rm2_line],
+    handles=[rm_line],
     fontsize=fs-1,
     loc="lower right",
     frameon=True,
     framealpha=0.85,
     facecolor="white"
 )
-
 plt.gca().add_artist(leg1)
 
-plt.grid(True)
+plt.grid(True, which='both', linestyle='--', linewidth=0.5)
 plt.xticks(fontsize=fs)
 plt.yticks(fontsize=fs)
 plt.tight_layout()
-plt.savefig("plots/BGO_transverse_profile_normed.pdf")
+plt.savefig("plots/BGO_transverse_profile_normed_log.pdf")
 
 # ============================================================
 # Plot 3: Transverse profile integrated over entire BGO
@@ -159,15 +159,14 @@ for E, lab, col in zip(energies, labels, colors):
 
 # 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 in cm", fontsize=fs)
 plt.ylabel(r"$\int \langle dE/(\mathrm{d}z\,\mathrm{d}A) \rangle\,dz$" + "\nin MeV/cm$^2$",fontsize=fs)
 plt.title("Transverse energy deposition integrated over 10 cm BGO (theory)", fontsize=fs+2)
 
-
+plt.xlim(0, 2.5)
-plt.xlim(0, r.max())
+plt.yscale("log")
-plt.ylim(0, None)
+plt.ylim(0.1, 100000)
 
 # --- Legend 1: Energies ---
 leg1 = plt.legend(
@@ -183,20 +182,19 @@ leg1 = plt.legend(
 
 # --- Legend 2: Geometry ---
 leg2 = plt.legend(
-    handles=[rm_line, rm2_line],
+    handles=[rm_line],
     fontsize=fs-1,
     loc="lower right",
     frameon=True,
     framealpha=0.85,
     facecolor="white"
 )
-
 plt.gca().add_artist(leg1)
 
-plt.grid(True)
+plt.grid(True, which='both', linestyle='--', linewidth=0.5)
 plt.xticks(fontsize=fs)
 plt.yticks(fontsize=fs)
 plt.tight_layout()
-plt.savefig("plots/BGO_transverse_profile_depth.pdf")
+plt.savefig("plots/BGO_transverse_profile_depth_log.pdf")
 
 plt.show()

shower_trans.py

@@ -28,7 +28,7 @@ colors = ["C0", "C1", "C2", "C3", "C4"]
 # -------------------------------------------------
 # Parameters
 # -------------------------------------------------
-r_max_cm = 8.0
+r_max_cm = 2.5
 n_bins   = 100
 fs = 26
 
@@ -84,54 +84,3 @@ plt.legend(title="Initial Energy", fontsize=fs-4, title_fontsize=fs-2, frameon=T
 plt.tight_layout()
 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_comp_slices.py

@@ -0,0 +1,180 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Transverse shower profiles in BGO
+Comparison: Theory (NKG, dashed) vs Geant4 simulation (step)
+
+Plot:
+- 5 slices (0-2,2-4,4-6,6-8,8-10 cm)
+- 1 sum plot (0-10 cm)
+- Legend: Initial energy + Theory vs Geant4
+- Slice-wise Integration für Theorie
+"""
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from scipy.stats import gamma
+import scienceplots
+from matplotlib.lines import Line2D
+
+# ============================================================
+# Style & Font
+# ============================================================
+plt.style.use(['science'])
+fs = 26
+
+# ============================================================
+# Material (BGO)
+# ============================================================
+X0 = 1.118
+Ec = 10.5
+b  = 0.5
+Rm = 2.26
+
+# ============================================================
+# Inputs
+# ============================================================
+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]
+labels   = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
+colors   = ["C0", "C1", "C2", "C3", "C4"]
+
+# ============================================================
+# Theory functions
+# ============================================================
+def tmax(E):
+    return np.log(E / Ec) - 0.5
+
+def dEdz(z, E):
+    t = z / X0
+    a = 1 + b * tmax(E)
+    return E * gamma.pdf(t, a, scale=1/b)  # MeV/cm
+
+# NKG transverse profile
+s = 1.0
+def rho_r(r):
+    x = r / Rm
+    core = np.zeros_like(r)
+    mask = r > 0
+    core[mask] = (x[mask]**(s-2)) * ((1 + x[mask])**(s-4.5))
+    if r[0] == 0:
+        core[0] = ((1e-9 / Rm)**(s-2)) * ((1 + 1e-9 / Rm)**(s-4.5))
+    norm = np.trapz(2*np.pi*r*core, r)
+    return core / norm
+
+# ============================================================
+# Radial binning
+# ============================================================
+r_max_cm = 2.5
+n_bins   = 80
+r_edges  = np.linspace(0, r_max_cm, n_bins+1)
+r_centers = 0.5*(r_edges[:-1]+r_edges[1:])
+A_ring = np.pi * (r_edges[1:]**2 - r_edges[:-1]**2)
+
+# ============================================================
+# Slices
+# ============================================================
+z_slices = [(0,2),(2,4),(4,6),(6,8),(8,10)]
+slice_labels = ["0-2 cm","2-4 cm","4-6 cm","6-8 cm","8-10 cm"]
+
+# ============================================================
+# Figure
+# ============================================================
+fig, axes = plt.subplots(6,1, figsize=(13,16), sharex=True,
+                         gridspec_kw={"height_ratios":[1,1,1,1,1,1],"hspace":0.05})
+
+z_grid = np.linspace(0,10,400)
+dz = z_grid[1]-z_grid[0]
+rho = rho_r(r_centers)
+
+# Handle-Listen für Energy-Legende (oben)
+energy_lines = []
+
+for file, E, lab, col in zip(files, energies, labels, colors):
+    df = pd.read_csv(file, sep="\t")
+    df = df[df["r"] <= r_max_cm]
+    N_events = df["event"].nunique()
+
+    # Für Startenergie-Legende
+    energy_lines.append(Line2D([0],[0], color=col, lw=2))
+
+    # --- Slices ---
+    for i, (z1, z2) in enumerate(z_slices):
+        ax = axes[i]
+
+        # Simulation
+        df_slice = df[(df["z"]>=z1)&(df["z"]<z2)]
+        E_sum_slice,_ = np.histogram(df_slice["r"], bins=r_edges, weights=df_slice["edep"])
+        sim_profile = E_sum_slice / (N_events*A_ring)
+
+        # Theorie Slice-wise
+        mask = (z_grid >= z1) & (z_grid < z2)
+        theory_profile = np.sum(dEdz(z_grid[mask], E)[:, None] * rho[None, :] * dz, axis=0)
+
+        ax.plot(r_centers, theory_profile, '--', lw=2, color=col)
+        ax.step(r_centers, sim_profile, where='mid', lw=2, color=col)
+
+        ax.set_yscale('log')
+        ax.set_ylim(0.05, 30000)
+        ax.set_ylabel(slice_labels[i], fontsize=fs-6)
+        ax.grid(True, which='both', ls='--', lw=0.6)
+
+    # --- Summe 0-10 cm ---
+    ax = axes[5]
+    # Simulation
+    E_sum, _ = np.histogram(df["r"], bins=r_edges, weights=df["edep"])
+    sim_profile = E_sum / (N_events * A_ring)
+    # Theorie über Gesamt 0-10 cm
+    theory_profile = np.sum(dEdz(z_grid, E)[:, None] * rho[None, :] * dz, axis=0)
+
+    ax.plot(r_centers, theory_profile, '--', lw=2, color=col)
+    ax.step(r_centers, sim_profile, where='mid', lw=2, color=col)
+
+    ax.set_yscale('log')
+    ax.set_ylim(0.05, 30000)
+    ax.set_ylabel('Sum', fontsize=fs-6)
+    ax.grid(True, which='both', ls='--', lw=0.6)
+
+# ============================================================
+# Cosmetics
+# ============================================================
+axes[-1].set_xlabel('Radius r in cm', fontsize=fs)
+for ax in axes:
+    ax.set_xlim(0,r_max_cm)
+    ax.tick_params(labelsize=fs-2)
+
+fig.text(0.04, 0.5, 
+         r"$\left\langle \int \frac{\mathrm{d}E}{\mathrm{d}z\,\mathrm{d}A} \, \mathrm{d}z \right\rangle$ in MeV/cm$^2$",
+         va='center', ha='center', rotation='vertical', fontsize=fs)
+
+# ============================================================
+# Überschrift
+# ============================================================
+fig.suptitle("Transverse shower profiles in BGO", fontsize=fs+2, y=0.94)
+
+# ============================================================
+# Legend: Initialenergie zwischen Titel und oberstem Plot
+# ============================================================
+fig.legend(energy_lines, labels, loc='upper center', bbox_to_anchor=(0.5, 0.92),
+           ncol=5, fontsize=fs-10, frameon=True)
+
+# ============================================================
+# Legend: Theory vs Geant4 Simulation
+# ============================================================
+sim_lines = [
+    Line2D([0],[0], color='k', lw=2, linestyle='-', label='Geant4 simulation'),
+    Line2D([0],[0], color='k', lw=2, linestyle='--', label='Theory')
+]
+axes[0].legend(handles=sim_lines, loc='upper right', ncol=2, fontsize=fs-10, frameon=True)
+
+plt.tight_layout(rect=[0,0,1,0.85])
+plt.savefig("plots/shower_trans_comp_slices.pdf", dpi=300)
+plt.show()

shower_trans_comparison.py

@@ -0,0 +1,224 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Transverse shower profiles in BGO
+Comparison: Theory (NKG, dashed) vs Geant4 simulation (step)
+Residuals shown logarithmically
+"""
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from scipy.stats import gamma
+from matplotlib.lines import Line2D
+import sys
+from pathlib import Path
+
+sys.path.append(str(Path(__file__).resolve().parents[1]))
+import plotstyle
+
+import scienceplots 
+
+# ============================================================
+# Style & Font
+# ============================================================
+plotstyle.set_style("paper")  # nutzt Slides-Mode
+fs = plt.rcParams["font.size"]
+
+# ============================================================
+# Material (BGO)
+# ============================================================
+X0 = 1.118
+Ec = 10.5
+b  = 0.5
+Rm = 2.26
+
+# ============================================================
+# Inputs
+# ============================================================
+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]
+labels   = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
+colors   = ["C0", "C1", "C2", "C3", "C4"]
+
+# ============================================================
+# Theory functions
+# ============================================================
+def tmax(E):
+    return np.log(E / Ec) - 0.5
+
+def dEdz(z, E):
+    t = z / X0
+    a = 1 + b * tmax(E)
+    return E * gamma.pdf(t, a, scale=1/b) / X0  # MeV/cm
+
+# NKG transverse profile
+s = 1.0
+def rho_r(r):
+    x = r / Rm
+    core = np.zeros_like(r)
+    mask = r > 0
+    core[mask] = (x[mask]**(s-2)) * ((1 + x[mask])**(s-4.5))
+    if r[0] == 0:
+        core[0] = ((1e-9 / Rm)**(s-2)) * ((1 + 1e-9 / Rm)**(s-4.5))
+    norm = np.trapz(2*np.pi*r*core, r)
+    return core / norm
+
+# ============================================================
+# Common binning
+# ============================================================
+r_max_cm = 2.5
+n_bins   = 80
+
+r_edges   = np.linspace(0, r_max_cm, n_bins + 1)
+r_centers = 0.5 * (r_edges[:-1] + r_edges[1:])
+r_in  = r_edges[:-1]
+r_out = r_edges[1:]
+A_ring = np.pi * (r_out**2 - r_in**2)
+
+# ============================================================
+# Legend styles
+# ============================================================
+style_legend = [
+    Line2D([0], [0], color="k", lw=2, linestyle="--", label="Theory"),
+    Line2D([0], [0], color="k", lw=2, linestyle="-", label="Geant4 simulation"),
+]
+
+# ============================================================
+# FIGURE 1 — Shower Maximum
+# ============================================================
+fig, (ax, ax_res) = plt.subplots(
+    2, 1,
+    sharex=True,
+    gridspec_kw={"height_ratios": [3, 1], "hspace": 0.08}
+)
+
+energy_lines = []
+dz_slice = 0.5
+
+for file, E, lab, col in zip(files, energies, labels, colors):
+
+    df = pd.read_csv(file, sep="\t")
+    N_events = df["event"].nunique()
+
+    z0 = tmax(E) * X0
+
+    df_slice = df[
+        (df["z"] >= z0 - dz_slice/2) &
+        (df["z"] <= z0 + dz_slice/2) &
+        (df["r"] <= r_max_cm)
+    ]
+
+    E_sum, _ = np.histogram(df_slice["r"], bins=r_edges, weights=df_slice["edep"])
+    sim_profile = E_sum / (N_events * dz_slice * A_ring)  # MeV/cm^2
+
+    # Theory
+    theory_profile = dEdz(z0, E) * rho_r(r_centers)  # MeV/cm^2
+
+    ax.plot(r_centers, theory_profile, linestyle="--", linewidth=2, color=col)
+    ax.step(r_centers, sim_profile, where="mid", linewidth=2, color=col)
+
+    energy_lines.append(Line2D([0], [0], color=col, lw=2, linestyle="-"))
+
+    # Residuals
+    ax_res.step(r_centers, (sim_profile - theory_profile) / theory_profile, color=col, linewidth=1.5)
+
+ax.set_yscale("log")
+ax.set_xlim(0, r_max_cm)
+ax.set_ylabel(
+    r"$\left\langle \frac{\mathrm{d}E}{\mathrm{d}z\,\mathrm{d}A} \right\rangle$ in MeV/cm$^2$",
+    fontsize=fs
+)
+ax.set_title("Transverse shower profile at shower maximum (BGO)")
+ax_res.set_xlim(0, r_max_cm)
+ax_res.set_xlabel("Radius r in cm")
+ax_res.set_ylabel(r"$\frac{\mathrm{G4sim} - \mathrm{Theory}}{\mathrm{Theory}}$",)
+ax_res.axhline(0, color='k', linestyle=':', linewidth=1.5)
+
+# Legends
+leg1 = ax.legend(
+    energy_lines, labels,
+    title="Initial Energy",
+    loc="upper right"
+)
+
+leg2 = ax.legend(
+    handles=style_legend,
+    loc="lower left",
+)
+ax.add_artist(leg1)
+
+plotstyle.savefig("shower_trans_comparison_max", category="BGO")
+
+# ============================================================
+# FIGURE 2 — Integrated over crystal
+# ============================================================
+fig, (ax, ax_res) = plt.subplots(
+    2, 1,
+    sharex=True,
+    gridspec_kw={"height_ratios": [3, 1], "hspace": 0.05}
+)
+
+z_grid = np.linspace(0, 10, 400)
+dz = z_grid[1] - z_grid[0]
+
+for file, E, lab, col in zip(files, energies, labels, colors):
+
+    df = pd.read_csv(file, sep="\t")
+    df = df[df["r"] <= r_max_cm]
+    N_events = df["event"].nunique()
+
+    # Simulation
+    E_sum, _ = np.histogram(df["r"], bins=r_edges, weights=df["edep"])
+    sim_profile = E_sum / (N_events * A_ring)  # MeV/cm^2
+
+    # Theory
+    prof_z = dEdz(z_grid, E)
+    theory_profile = np.sum(
+        prof_z[:, None] * rho_r(r_centers)[None, :] * dz,
+        axis=0
+    )
+
+    ax.plot(r_centers, theory_profile, linestyle="--", linewidth=2, color=col)
+    ax.step(r_centers, sim_profile, where="mid", linewidth=2, color=col)
+    ax_res.step(r_centers, (sim_profile - theory_profile)/theory_profile, color=col, linewidth=1.5)
+
+ax.set_yscale("log")
+ax.set_ylim(0.05, 30000)
+ax.set_xlim(0, r_max_cm)
+ax.set_ylabel(
+    r"$\left\langle \int \frac{\mathrm{d}E}{\mathrm{d}z\,\mathrm{d}A} \, \mathrm{d}z \right\rangle$ in MeV/cm$^2$",
+)
+ax.set_title("Transverse shower profile integrated over 10 cm BGO")
+
+ax_res.set_xlim(0, r_max_cm)
+ax_res.set_xlabel("Radius r in cm")
+ax_res.set_ylabel(
+    r"$\frac{\mathrm{G4sim} - \mathrm{Theory}}{\mathrm{Theory}}$",
+)
+ax_res.axhline(0, color='k', linestyle=':', linewidth=1.5)
+
+# Legends
+leg1 = ax.legend(
+    energy_lines, labels,
+    title="Initial Energy",
+    title_fontsize=fs-2,
+    loc="upper right",
+)
+
+leg2 = ax.legend(
+    handles=style_legend,
+    loc="lower left",
+)
+ax.add_artist(leg1)
+
+plotstyle.savefig("shower_trans_comparison_integrated", category="BGO")
+
+plt.show()

shower_trans_max.py

@@ -39,7 +39,7 @@ colors   = ["C0", "C1", "C2", "C3", "C4"]
 # -------------------------------------------------
 # Analysis parameters
 # -------------------------------------------------
-r_max_cm = 8.0
+r_max_cm = 2.5
 n_bins   = 80
 
 dz_slice = 0.5   # cm  (slice thickness around shower maximum)
@@ -111,11 +111,113 @@ 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.grid(True)
 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.tight_layout()
 plt.savefig("plots/G4_transverse_profile_showermax.png", dpi=300)
+
+# --------------------------
+# Fixed z positions
+# --------------------------
+z_positions = [2.0, 4.0, 6.0]  # cm
+dz_slice = 0.5                # cm slice thickness
+n_slices = len(z_positions)
+
+# --------------------------
+# Figure
+# --------------------------
+plt.style.use(['science'])
+fig, axes = plt.subplots(n_slices, 1, figsize=(10, 10), sharex=True)
+plt.subplots_adjust(hspace=0.12)
+
+fig.suptitle(
+    "Transverse energy density at fixed depths in BGO (Geant4)",
+    fontsize=fs+2,
+    y=0.95
+)
+
+# Dummy lines for legend (energy only)
+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,
+    loc="upper center",
+    bbox_to_anchor=(0.5, 1.25)
+)
+
+# Shared y-label
+fig.text(
+    0.02, 0.5,
+    r"$\langle \mathrm{d}E/(\mathrm{d}z\,\mathrm{d}A) \rangle$ in MeV/cm$^3$",
+    va="center",
+    rotation="vertical",
+    fontsize=fs
+)
+
+# --------------------------
+# Loop over z slices
+# --------------------------
+for ax, z_ref in zip(axes, z_positions):
+
+    for file, E, color in zip(files, energies, colors):
+        df = pd.read_csv(file, sep="\t")
+
+        # Number of primary events
+        N_events = df["event"].nunique()
+
+        # z slice around fixed depth
+        df_slice = df[
+            (df["z"] >= z_ref - dz_slice/2) &
+            (df["z"] <= z_ref + dz_slice/2)
+        ]
+
+        df_slice = df_slice[df_slice["r"] <= r_max_cm]
+
+        # Radial histogram
+        r_edges = np.linspace(0.0, r_max_cm, n_bins + 1)
+        r_centers = 0.5 * (r_edges[:-1] + r_edges[1:])
+
+        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)
+
+        # Normalised profile
+        profile = E_sum / (N_events * dz_slice * A_ring)
+
+        ax.step(
+            r_centers,
+            profile,
+            where="mid",
+            lw=2,
+            color=color
+        )
+
+    ax.set_xlim(0, r_max_cm)
+    ax.set_yscale("log")
+    ax.grid(True, ls="--", lw=0.6)
+    ax.tick_params(labelsize=fs-2)
+
+    # z-label on the left of each panel
+    ax.set_ylabel(f"z = {z_ref:.0f} cm", fontsize=fs)
+
+# Bottom x-label
+axes[-1].set_xlabel("Radius r in cm", fontsize=fs)
+
+plt.tight_layout(rect=[0.06, 0.03, 0.98, 0.93])
+plt.savefig("plots/G4_transverse_profile_fixed_z.png", dpi=300)
 plt.show()
+

shower_trans_slices.py

@@ -9,6 +9,10 @@ import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
 import scienceplots 
+from matplotlib.lines import Line2D
+
+plt.style.use(['science'])
+fs = 26
 
 # --------------------------
 # Input files and parameters
@@ -30,7 +34,7 @@ fs = 24
 # Analysis parameters
 # --------------------------
 z_max_cm = 10.0
-r_max_cm = 8.0
+r_max_cm = 2.5
 dz_layer  = 2.0   # cm
 n_bins_r  = 100
 
@@ -44,20 +48,8 @@ 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)
-
-# 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,
-               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$",
-         va='center', rotation='vertical', fontsize=fs)
 
 # --------------------------
 # Layerwise plotting with integrals
@@ -96,8 +88,8 @@ 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.01, 1e5)
-    ax.grid(True)
+    ax.grid(True, which='both', linestyle='--', linewidth=0.5)
     ax.tick_params(labelsize=fs-2)
 
     if i < n_layers:
@@ -109,83 +101,26 @@ for i, (z_min, z_max) in enumerate(layers + [(0,z_max_cm)]):
     # Legende innerhalb der Achse
     ax.legend(fontsize=12, frameon=True, framealpha=0.85, facecolor="white")
 
+# Cosmetics
+axes[-1].set_xlabel('Radius r in cm', fontsize=fs)
+for ax in axes:
+    ax.set_xlim(0,r_max_cm)
+    ax.tick_params(labelsize=fs-2)
+
+# Gemeinsames y-Label
+fig.text(0.04, 0.5,
+         r"$\left\langle \int \frac{dE}{dz\,dA} dz \right\rangle$ [MeV/cm$^2$]",
+         va='center', ha='center', rotation='vertical', fontsize=fs)
+
+# Überschrift
+fig.suptitle("Transverse shower profiles in BGO (Geant4 Simulation)", fontsize=fs+2, y=0.98)
+
+# Startenergie-Legende oberhalb oberstem Plot
+energy_lines = [Line2D([0],[0], color=c,lw=2) for c in colors]
+fig.legend(energy_lines, labels, loc='upper center', bbox_to_anchor=(0.5,0.95),
+           ncol=5, fontsize=fs-10, frameon=True)
+
 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()
+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_slices2.py

@@ -0,0 +1,121 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Transverse shower profiles in BGO
+Nur Geant4 Simulation
+Legenden: obere Legende = Startenergien, rechts pro Subplot = integrierte Energie pro Slice
+"""
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import scienceplots
+from matplotlib.lines import Line2D
+
+plt.style.use(['science'])
+fs = 26
+
+# Radial binning
+r_max_cm = 2.5
+n_bins = 80
+r_edges = np.linspace(0, r_max_cm, n_bins+1)
+r_centers = 0.5*(r_edges[:-1] + r_edges[1:])
+r_in = r_edges[:-1]
+r_out = r_edges[1:]
+A_ring = np.pi*(r_out**2 - r_in**2)
+
+# Slices
+z_slices = [(0,2), (2,4), (4,6), (6,8), (8,10)]
+slice_labels = ["0-2 cm", "2-4 cm", "4-6 cm", "6-8 cm", "8-10 cm"]
+
+# Inputs
+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"]
+
+# Figure
+fig, axes = plt.subplots(6, 1, figsize=(13,16), sharex=True,
+                         gridspec_kw={"height_ratios":[1,1,1,1,1,1], "hspace":0.05})
+
+# Bereite Listen für Handles und Labels pro Subplot vor (nur integrierte Energie)
+handles_slices = [[] for _ in range(len(z_slices)+1)]
+labels_slices = [[] for _ in range(len(z_slices)+1)]
+
+# Handle-Liste für Startenergie-Legende oben
+energy_lines = []
+
+for file, E, col, label in zip(files, [100,200,500,1000,2000], colors, labels):
+    df = pd.read_csv(file, sep="\t")
+    df = df[df["r"] <= r_max_cm]
+    N_events = df["event"].nunique()
+
+    # Für Startenergie-Legende
+    energy_lines.append(Line2D([0], [0], color=col, lw=2))
+
+    for i, (z1, z2) in enumerate(z_slices):
+        ax = axes[i]
+        df_slice = df[(df["z"] >= z1) & (df["z"] < z2)]
+        E_sum_slice, _ = np.histogram(df_slice["r"], bins=r_edges, weights=df_slice["edep"])
+        sim_profile = E_sum_slice / (N_events * A_ring)
+
+        # Physikalisch korrekte Slice-Energie
+        E_slice_total = np.sum(E_sum_slice) / N_events
+        E_slice_total = min(E_slice_total, E)
+
+        line = ax.step(r_centers, sim_profile, where='mid', lw=2, color=col)[0]
+
+        # Nur integrierte Energie als Subplot-Legende
+        handles_slices[i].append(line)
+        labels_slices[i].append(f"{E_slice_total:.1f} MeV")
+
+        ax.set_yscale('log')
+        ax.set_ylim(0.05, 30000)
+        ax.set_ylabel(slice_labels[i], fontsize=fs-6)
+        ax.grid(True, which='both', ls='--', lw=0.6)
+
+    # Summe 0-10 cm
+    ax = axes[5]
+    E_sum, _ = np.histogram(df["r"], bins=r_edges, weights=df["edep"])
+    sim_profile = E_sum / (N_events * A_ring)
+    E_total = np.sum(E_sum) / N_events
+    E_total = min(E_total, E)
+
+    line = ax.step(r_centers, sim_profile, where='mid', lw=2, color=col)[0]
+    handles_slices[5].append(line)
+    labels_slices[5].append(f"{E_total:.1f} MeV")
+
+    ax.set_yscale('log')
+    ax.set_ylim(0.05, 30000)
+    ax.set_ylabel('Sum', fontsize=fs-6)
+    ax.grid(True, which='both', ls='--', lw=0.6)
+
+# Füge Subplot-Legenden (nur integrierte Energie) hinzu
+for ax, h, l in zip(axes, handles_slices, labels_slices):
+    ax.legend(h, l, loc='upper right', fontsize=fs-9, frameon=True)
+
+# Cosmetics
+axes[-1].set_xlabel('Radius r in cm', fontsize=fs)
+for ax in axes:
+    ax.set_xlim(0, r_max_cm)
+    ax.tick_params(labelsize=fs-2)
+
+# Gemeinsames y-Label
+fig.text(0.04, 0.5, r"$\left\langle \int \frac{dE}{dz\,dA} dz \right\rangle$ in MeV/cm$^2$",
+         va='center', ha='center', rotation='vertical', fontsize=fs)
+
+# Überschrift
+fig.suptitle("Transverse shower profiles in BGO (Geant4 Simulation)", fontsize=fs+2, y=0.94)
+
+# Startenergie-Legende oben zwischen Titel und oberstem Plot
+fig.legend(energy_lines, labels, loc='upper center', bbox_to_anchor=(0.5,0.92),
+           ncol=5, fontsize=fs-10, frameon=True)
+
+plt.tight_layout(rect=[0,0,1,0.85])
+plt.savefig("plots/G4_transverse_layers_sum2.pdf", dpi=300)
+plt.show()

shower_trans_theory.py

@@ -1,13 +1,17 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
-transversal shower profile theory slicewise
+transversal shower profile theory slicewise (BGO, NKG-corrected)
 """
 
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.stats import gamma
 import scienceplots 
+from matplotlib.lines import Line2D
+
+plt.style.use(['science'])
+fs = 26
 
 # --------------------------
 # Material & Parameter
@@ -20,8 +24,6 @@ Rm = 2.26
 energies = [100, 200, 500, 1000, 2000]  # MeV
 colors = ["C0","C1","C2","C3","C4"]
 
-fs = 24  # globale Schriftgröße
-
 # --------------------------
 # Profile
 # --------------------------
@@ -30,15 +32,41 @@ def dEdz(z, E):
     a = 1 + b*(np.log(E/Ec)-0.5)
     return E * gamma.pdf(t, a, scale=1/b) / X0  # MeV/cm
 
-def rho_r(r):
+# --------------------------
-    return np.exp(-r**2/(2*Rm**2)) / (2*np.pi*Rm**2)
+# NKG transverse profile (safe at r=0)
+# --------------------------
+s = 1.0
+
+# Verwende r_norm für Normierung (0–8 cm)
+r_norm = np.linspace(0, 2.5, 300)
+
+def rho_r(r, r_norm=r_norm):
+    """
+    NKG transverse profile
+    normalized: ∫ 2π r ρ(r) dr = 1
+    avoids divide-by-zero at r=0
+    """
+    x = r / Rm
+    core = np.zeros_like(r)
+    mask = r > 0
+    core[mask] = (x[mask]**(s-2)) * ((1 + x[mask])**(s-4.5))
+    if r[0] == 0:
+        core[0] = ((1e-9 / Rm)**(s-2)) * ((1 + 1e-9 / Rm)**(s-4.5))
+    # Normiere über r_norm (0–8 cm)
+    x_norm = r_norm / Rm
+    core_norm = np.zeros_like(r_norm)
+    mask_norm = r_norm > 0
+    core_norm[mask_norm] = (x_norm[mask_norm]**(s-2)) * ((1 + x_norm[mask_norm])**(s-4.5))
+    core_norm[0] = ((1e-9 / Rm)**(s-2)) * ((1 + 1e-9 / Rm)**(s-4.5))
+    norm = np.trapz(2*np.pi*r_norm*core_norm, r_norm)
+    return core / norm
 
 # --------------------------
 # Gitter
 # --------------------------
 z_grid = np.linspace(0, 10, 400)
 dz = z_grid[1]-z_grid[0]
-r_grid = np.linspace(0, 8, 300)
+r_grid = np.linspace(0, 2.5, 300)  # nur bis 2.5 cm
 dr = r_grid[1]-r_grid[0]
 
 # --------------------------
@@ -50,24 +78,13 @@ 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)
-
-# 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,
-               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)
 
 # --------------------------
 # Plots
 # --------------------------
+
 for i, (z_min, z_max) in enumerate(layers + [(0,10)]):
     ax = axes[i]
     z_mask = (z_grid >= z_min) & (z_grid < z_max)
@@ -80,9 +97,9 @@ for i, (z_min, z_max) in enumerate(layers + [(0,10)]):
         linestyle = '--' if i==n_layers else '-'  # Summe gestrichelt
         ax.plot(r_grid, radial_profile, color=color, linewidth=2, linestyle=linestyle, label=f"{E_layer:.1f} MeV")
     
-    ax.set_xlim(0, 8)
+    ax.set_xlim(0, 2.5)
-    ax.set_ylim(0, None)
+    ax.set_ylim(0.01, 1e5)
-    ax.grid(True)
+    ax.grid(True, which='both', linestyle='--', linewidth=0.5)
     ax.tick_params(labelsize=fs-2)
     
     if i < n_layers:
@@ -91,8 +108,28 @@ for i, (z_min, z_max) in enumerate(layers + [(0,10)]):
         ax.set_ylabel("Sum", fontsize=fs)
         ax.set_xlabel("Radius r in cm", fontsize=fs)
     
+    ax.set_yscale("log")  # logarithmische y-Achse
     ax.legend(fontsize=12, frameon=True, framealpha=0.85, facecolor="white")
 
+# Cosmetics
+axes[-1].set_xlabel('Radius r in cm', fontsize=fs)
+for ax in axes:
+    ax.set_xlim(0,r_max_cm)
+    ax.tick_params(labelsize=fs-2)
+
+# Gemeinsames y-Label
+fig.text(0.04, 0.5,
+         r"$\left\langle \int \frac{dE}{dz\,dA} dz \right\rangle$ [MeV/cm$^2$]",
+         va='center', ha='center', rotation='vertical', fontsize=fs)
+
+# Überschrift
+fig.suptitle("Transverse shower profiles in BGO (Theory)", fontsize=fs+2, y=0.98)
+
+# Startenergie-Legende oberhalb oberstem Plot
+energy_lines = [Line2D([0],[0], color=c,lw=2) for c in colors]
+fig.legend(energy_lines, labels, loc='upper center', bbox_to_anchor=(0.5,0.95),
+           ncol=5, fontsize=fs-10, frameon=True)
+
 plt.tight_layout(rect=[0.05,0.03,0.97,0.93])
 plt.savefig("plots/shower_transverse_layers_sum.pdf")
 plt.savefig("plots/shower_transverse_layers_sum.png", dpi=300)

shower_trans_theory2.py

@@ -0,0 +1,133 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Transverse shower profiles in BGO
+Nur Theorie
+Legenden: obere Legende = Startenergien, rechts pro Subplot = integrierte Energie pro Slice
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.stats import gamma
+import scienceplots
+from matplotlib.lines import Line2D
+
+plt.style.use(['science'])
+fs = 26
+
+# Material
+X0 = 1.118
+Ec = 10.5
+b  = 0.5
+Rm = 2.26
+
+# Energies & Farben
+energies = [100, 200, 500, 1000, 2000]
+labels = ["100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
+colors = ["C0", "C1", "C2", "C3", "C4"]
+
+# Theorie-Funktionen
+def tmax(E):
+    return np.log(E/Ec)-0.5
+
+def dEdz(z, E):
+    t = z/X0
+    a = 1 + b*tmax(E)
+    return E*gamma.pdf(t, a, scale=1/b)  # MeV/cm
+
+def rho_r(r):
+    s = 1.0
+    x = r/Rm
+    core = np.zeros_like(r)
+    mask = r > 0
+    core[mask] = (x[mask]**(s-2)) * ((1+x[mask])**(s-4.5))
+    if r[0] == 0:
+        core[0] = ((1e-9/Rm)**(s-2)) * ((1+1e-9/Rm)**(s-4.5))
+    norm = np.trapz(2*np.pi*r*core, r)
+    return core/norm
+
+# Radial & Slices
+r_max_cm = 2.5
+n_bins = 80
+r_edges = np.linspace(0, r_max_cm, n_bins+1)
+r_centers = 0.5*(r_edges[:-1] + r_edges[1:])
+r_in = r_edges[:-1]
+r_out = r_edges[1:]
+A_ring = np.pi*(r_out**2 - r_in**2)
+
+z_slices = [(0,2), (2,4), (4,6), (6,8), (8,10)]
+slice_labels = ["0-2 cm", "2-4 cm", "4-6 cm", "6-8 cm", "8-10 cm"]
+
+z_grid = np.linspace(0, 10, 400)
+dz = z_grid[1] - z_grid[0]
+rho = rho_r(r_centers)
+
+# Figure
+fig, axes = plt.subplots(6, 1, figsize=(13,16), sharex=True,
+                         gridspec_kw={"height_ratios":[1,1,1,1,1,1], "hspace":0.05})
+
+# Listen für Subplot-Legenden vorbereiten
+handles_slices = [[] for _ in range(len(z_slices)+1)]
+labels_slices = [[] for _ in range(len(z_slices)+1)]
+
+for E, lab, col in zip(energies, labels, colors):
+    for i, (z1, z2) in enumerate(z_slices):
+        ax = axes[i]
+        mask = (z_grid >= z1) & (z_grid < z2)
+        theory_profile = np.sum(dEdz(z_grid[mask], E)[:,None]*rho[None,:]*dz, axis=0)
+
+        # Slice-Energie
+        E_slice_total = np.sum(theory_profile*A_ring)
+        E_slice_total = min(E_slice_total, E)
+
+        line, = ax.plot(r_centers, theory_profile, lw=2, color=col)
+
+        handles_slices[i].append(line)
+        labels_slices[i].append(f"{E_slice_total:.1f} MeV")  # nur integrierte Energie
+
+        ax.set_yscale('log')
+        ax.set_ylim(0.05, 30000)
+        ax.set_ylabel(slice_labels[i], fontsize=fs-6)
+        ax.grid(True, which='both', ls='--', lw=0.6)
+
+    # Summe 0-10 cm
+    ax = axes[5]
+     # Theorie Slice-wise
+    mask = (z_grid >= z1) & (z_grid < z2)
+    theory_profile = np.sum(dEdz(z_grid[mask], E)[:, None] * rho[None, :] * dz, axis=0)
+    E_total = np.sum(theory_profile*A_ring)
+    E_total = min(E_total, E)
+
+    line, = ax.plot(r_centers, theory_profile, lw=2, color=col)
+    handles_slices[5].append(line)
+    labels_slices[5].append(f"{E_total:.1f} MeV")
+
+    ax.set_yscale('log')
+    ax.set_ylim(0.05, 30000)
+    ax.set_ylabel('Sum', fontsize=fs-6)
+    ax.grid(True, which='both', ls='--', lw=0.6)
+
+# Subplot-Legenden (integrierte Energie pro Slice)
+for ax, h, l in zip(axes, handles_slices, labels_slices):
+    ax.legend(h, l, loc='upper right', fontsize=fs-9, frameon=True)
+
+# Cosmetics
+axes[-1].set_xlabel('Radius r in cm', fontsize=fs)
+for ax in axes:
+    ax.set_xlim(0, r_max_cm)
+    ax.tick_params(labelsize=fs-2)
+
+fig.text(0.04, 0.5,
+         r"$\left\langle \int \frac{dE}{dz\,dA} dz \right\rangle$ in MeV/cm$^2$",
+         va='center', ha='center', rotation='vertical', fontsize=fs)
+
+fig.suptitle("Transverse shower profiles in BGO (Theory)", fontsize=fs+2, y=0.94)
+
+# Startenergie-Legende oben zwischen Titel und oberstem Plot
+energy_lines = [Line2D([0],[0], color=c, lw=2) for c in colors]
+fig.legend(energy_lines, labels, loc='upper center', bbox_to_anchor=(0.5,0.92),
+           ncol=5, fontsize=fs-10, frameon=True)
+
+plt.tight_layout(rect=[0,0,1,0.85])
+plt.savefig("plots/shower_transverse_layers_sum2.pdf", dpi=300)
+plt.show()

showermap_comparison.py

@@ -0,0 +1,222 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+6-Panel plot:
+- Top-left: Theory 2D contours (absolute dE/dz in MeV/cm²)
+- Other 5: Geant4 heatmaps for 100,200,500,1GeV,2GeV
+- Shared colorbar for Geant4 heatmaps, vmax=0.8 MeV/cm²
+- Subcaptions: "Theory", "100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"
+- Theory contours overlayed on Geant4 heatmaps
+"""
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from matplotlib.gridspec import GridSpec
+from matplotlib.lines import Line2D
+from scipy.stats import gamma
+import scienceplots
+
+plt.style.use(['science'])
+fs = 22
+
+# --------------------------
+# Material & BGO parameters
+# --------------------------
+X0_cm = 1.118
+Ec_MeV = 10.5
+b = 0.5
+Rm_cm = 2.26
+max_depth_cm = 10
+
+# --------------------------
+# Energies
+# --------------------------
+energies_MeV = [100, 200, 500, 1000, 2000]
+colors = ["C0", "C1", "C2", "C3", "C4"]
+
+# --------------------------
+# Theory functions
+# --------------------------
+def tmax(E, Ec):
+    return np.log(E / Ec) - 0.5
+
+def dEdz(z, E, X0=X0_cm, b=0.5, Ec=Ec_MeV):
+    """Longitudinal deposition MeV/cm"""
+    t = z / X0
+    a = 1 + b * tmax(E, Ec)
+    return E * gamma.pdf(t, a, scale=1 / b)  # MeV/cm
+
+s = 1.0
+def rho_r(r, Rm=Rm_cm, s_val=1.0):
+    x = r / Rm
+    core = np.zeros_like(r)
+    mask = r > 0
+    core[mask] = (x[mask] ** (s_val - 2)) * ((1 + x[mask]) ** (s_val - 4.5))
+    if r[0] == 0:
+        core[0] = ((1e-9 / Rm) ** (s_val - 2)) * ((1 + 1e-9 / Rm) ** (s_val - 4.5))
+    norm = np.trapz(2 * np.pi * r * core, r)
+    return core / norm
+
+# --------------------------
+# Radial & depth grid
+# --------------------------
+r_max = 2.5
+r = np.linspace(0, r_max, 300)
+z = np.linspace(0, max_depth_cm, 400)
+Z, R = np.meshgrid(z, r, indexing="ij")
+
+# --------------------------
+# Geant4 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",
+]
+
+# --------------------------
+# Figure layout
+# --------------------------
+fig = plt.figure(figsize=(18, 10))
+gs = GridSpec(2, 3, figure=fig, wspace=0.3, hspace=0.4)
+
+axes = [
+    fig.add_subplot(gs[0, 0]),  # Theory contours
+    fig.add_subplot(gs[0, 1]),  # 100 MeV
+    fig.add_subplot(gs[0, 2]),  # 200 MeV
+    fig.add_subplot(gs[1, 0]),  # 500 MeV
+    fig.add_subplot(gs[1, 1]),  # 1 GeV
+    fig.add_subplot(gs[1, 2]),  # 2 GeV
+]
+
+subtitles = ["Theory", "100 MeV", "200 MeV", "500 MeV", "1 GeV", "2 GeV"]
+
+# --------------------------
+# 1) Theory contour plot
+# --------------------------
+ax = axes[0]
+heatmap_levels = [0.1, 1, 10]  # MeV/cm²
+linestyles = [":", "--", "-"]
+
+for E, col in zip(energies_MeV, colors):
+    long_prof = dEdz(z, E)
+    shower_2D = np.zeros((len(z), len(r)))
+    rho = rho_r(r)
+    rho[0] = rho[1]
+    for i in range(len(z)):
+        shower_2D[i, :] = long_prof[i] * rho
+
+    ax.contour(
+        R, Z, shower_2D,
+        levels=heatmap_levels,
+        colors=[col] * len(heatmap_levels),
+        linestyles=linestyles,
+        linewidths=1.5
+    )
+
+ax.set_xlabel("Radius r in cm", fontsize=fs)
+ax.set_ylabel("Depth z in cm", fontsize=fs)
+ax.set_xlim(0, r_max)
+ax.set_ylim(max_depth_cm, 0)
+ax.set_title("2D Energy density contours", fontsize=fs)
+ax.grid(True, linestyle="--", linewidth=0.5)
+ax.tick_params(labelsize=fs)
+
+line_legend = [Line2D([0], [0], color="k", ls=ls, lw=2) for ls in linestyles]
+ax.legend(
+    line_legend,
+    [f"{lvl} MeV/cm²" for lvl in heatmap_levels],
+    title="Contour level",
+    fontsize=fs - 2,
+    title_fontsize=fs,
+    loc="lower right",
+    frameon=True,
+)
+
+# --------------------------
+# 2–6) Geant4 heatmaps + Theory overlay
+# --------------------------
+vmin_cb = 0.0
+vmax_cb = 0.8
+im_list = []
+
+for ax, file, subtitle, col in zip(axes[1:], files, subtitles[1:], colors):
+    df = pd.read_csv(file, sep="\t")
+    df = df[df["r"] <= r_max]
+
+    dz = 0.1
+    dr = 0.1
+    z_edges = np.arange(0, max_depth_cm + dz, dz)
+    r_edges = np.arange(0, r_max + dr, dr)
+
+    H, _, _ = np.histogram2d(
+        df["z"], df["r"],
+        bins=[z_edges, r_edges],
+        weights=df["edep"]
+    )
+    counts, _, _ = np.histogram2d(
+        df["z"], df["r"],
+        bins=[z_edges, r_edges]
+    )
+
+    H_avg = np.zeros_like(H)
+    mask = counts > 0
+    H_avg[mask] = H[mask] / counts[mask]
+
+    Z_grid, R_grid = np.meshgrid(
+        (z_edges[:-1] + z_edges[1:]) / 2,
+        (r_edges[:-1] + r_edges[1:]) / 2,
+        indexing="ij"
+    )
+
+    im = ax.pcolormesh(
+        R_grid, Z_grid, H_avg,
+        shading="auto",
+        cmap="viridis",
+        vmin=vmin_cb,
+        vmax=vmax_cb
+    )
+    im_list.append(im)
+
+    # Theory overlay contours
+    E_val = float(subtitle.split()[0])
+    long_prof = dEdz(z, E_val)
+    shower_2D = np.zeros((len(z), len(r)))
+    rho = rho_r(r)
+    rho[0] = rho[1]
+    for i in range(len(z)):
+        shower_2D[i, :] = long_prof[i] * rho
+
+    ax.contour(
+        R, Z, shower_2D,
+        levels=heatmap_levels,
+        colors="white",
+        linestyles=linestyles,
+        linewidths=1
+    )
+
+    ax.set_title(subtitle, fontsize=fs, color=col, pad=10)
+    ax.set_xlabel("Radius r in cm", fontsize=fs)
+    ax.set_ylabel("Depth z in cm", fontsize=fs)
+    ax.set_ylim(max_depth_cm, 0)
+    ax.set_xlim(0, r_max)
+    ax.grid(True, linestyle="--", linewidth=0.5)
+    ax.tick_params(labelsize=fs)
+
+# --------------------------
+# Shared colorbar
+# --------------------------
+ymax = axes[1].get_position().y1
+ymin = axes[-1].get_position().y0
+cbar_ax = fig.add_axes([0.93, ymin, 0.015, ymax - ymin])
+cbar = fig.colorbar(im_list[0], cax=cbar_ax)
+cbar.set_label("dE/dz in MeV/cm²", fontsize=fs)
+cbar.ax.tick_params(labelsize=fs - 2)
+
+plt.tight_layout(rect=[0, 0, 0.92, 0.95])
+plt.savefig("plots/shower_map_comparison.pdf", dpi=300)
+plt.show()