xxx

BB-chaos.py

@@ -8,6 +8,7 @@ Created on Wed Feb 12 09:20:07 2025
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.constants import *
+import pandas as pd
 
 e0 = epsilon_0
 me = m_e
@@ -23,20 +24,20 @@ u = 1e-3 / Na
 
 particles = {
     'Proton': {
-        'E0_min': 50,   # MeV
+        'E0_min': 10,   # MeV
-        'E0_max': 2000, # MeV
+        'E0_max': 3000, # MeV
         'm0': m_p,
         'z': 1
     },
     'Helium': {
-        'E0_min': 1000,  # MeV
+        'E0_min': 100,  # MeV
-        'E0_max': 3000,  # MeV
+        'E0_max': 5000,  # MeV
         'm0': 6.644657e-27,
         'z': 2
     },
     'Muon': {
-        'E0_min': 50,    # MeV
+        'E0_min': 10,    # MeV
-        'E0_max': 2000,  # MeV
+        'E0_max': 3000,  # MeV
         'm0': 1.883531627e-28,
         'z': 1
     }
@@ -121,174 +122,221 @@ def bethebloch(Ekin, m0, z, ne, I):
     C = ne * z**2 * e**4 / (4 * np.pi * me * v**2 * e0**2)
     D = np.log(2 * lorentz_gamma**2 * me * v**2 / I)
     beta = v**2 / c**2
-    return -C * (D - beta)
+    return C * (D - beta)
 
-dx = 1e-9
+
+dx = 1e-7
 
 def Eloss(E0, m0, z, ne, I, dx, distance):
     steps = round(distance / dx)
     E = MeV_to_J(E0)
     for i in range(steps):
         dE = bethebloch(E, m0, z, ne, I) * dx
-        E += dE
+        E -= dE
     Ediff = MeV_to_J(E0) - E
-    return J_to_MeV(Ediff), J_to_MeV(E)
+    return J_to_MeV(Ediff)
 
-#Energieverluste berechnen und ausgeben für Si und BGO
+######################################################################################################################
-for material, params in materials.items():
+
-    print(f"--- Energieverluste in {material} ---")
+def calculate_Elosses(materials, particles):
-    for name, particle_params in particles.items():
+    for material, params in materials.items():
-        Ediff_min, E_min = Eloss(particle_params['E0_min'], particle_params['m0'], particle_params['z'], params['ne'], params['I'], dx, params['d'])
+        print(f"--- Energieverluste in {material} ---")
-        Ediff_max, E_max = Eloss(particle_params['E0_max'], particle_params['m0'], particle_params['z'], params['ne'], params['I'], dx, params['d'])
+        for name, particle_params in particles.items():
-        print(f'Eloss {name} in {material}:')
+            Ediff_min, E_min = Eloss(particle_params['E0_min'], particle_params['m0'], particle_params['z'], params['ne'], params['I'], dx, params['d'])
-        print(f"E0 = {particle_params['E0_min']} MeV      |  {particle_params['E0_max']} MeV")
+            Ediff_max, E_max = Eloss(particle_params['E0_max'], particle_params['m0'], particle_params['z'], params['ne'], params['I'], dx, params['d'])
-        print(f"E  = {E_min:.4f} MeV   |  E  = {E_max:.4f} MeV")
+            print(f'Eloss {name} in {material}:')
-        print(f"Ediff = {Ediff_min:.4f} MeV | Ediff = {Ediff_max:.4f} MeV")
+            print(f"E0 = {particle_params['E0_min']} MeV      |  {particle_params['E0_max']} MeV")
-        print('--------------------------')
+            print(f"E  = {E_min:.4f} MeV   |  E  = {E_max:.4f} MeV")
+            print(f"Ediff = {Ediff_min:.4f} MeV | Ediff = {Ediff_max:.4f} MeV")
+            print('--------------------------')
         
 ######################################################################################################################
 
-# Linienstile für alle Materialien (jetzt Schwarz für das Material)
-linestyles = {
-    'Si': '-',         # Linienstil für Si
-    'Bi': '--',        # Linienstil für Bi
-    'Ge': ':',         # Linienstil für Ge
-    'O': '-.',         # Linienstil für O
-    'BGO': '-'         # Linienstil für BGO (nun geändert)
-}
 
-# Farben für die Teilchen
+def plot_Eloss_E():
-colors = {
+    linestyles = {
-    'Proton': 'red',
+        'Si': '-',         # Linienstil für Si
-    'Helium': 'green',
+        'Bi': '--',        # Linienstil für Bi
-    'Muon': (0.5, 0.3, 0.8)  # RGB-Farbdefinition für Muon
+        'Ge': ':',         # Linienstil für Ge
-}
+        'O': '-.',         # Linienstil für O
-
+        'BGO': '-'         # Linienstil für BGO (nun geändert)
-# Erstelle den Plot für alle 3 Teilchen in verschiedenen Materialien
+    }
-plt.figure(figsize=(10,6))
+    
-
+    colors = {
-# Listen für die Legende
+        'Proton': 'red',
-legend_particle_handles = []
+        'Helium': 'green',
-legend_particle_labels = []
+        'Muon': (0.5, 0.3, 0.8)  # RGB-Farbdefinition für Muon
-legend_material_handles = []
+    }
-legend_material_labels = []
+    
-
+    plt.figure(figsize=(10,6))
-# dE/dx vs. Energie (in MeV) für alle Materialien
+    
-for material, params in materials.items():
+    legend_particle_handles = []
-    for name, particle_params in particles.items():
+    legend_particle_labels = []
-        E = np.linspace(particle_params['E0_min'], particle_params['E0_max'], 100000) * 1e6 * e  # Umrechnung in Joule
+    legend_material_handles = []
-        dEdx = -bethebloch(E, particle_params['m0'], particle_params['z'], params['ne'], params['I'])
+    legend_material_labels = []
+    
+    # dE/dx vs. Energie (in MeV) für alle Materialien
+    for material, params in materials.items():
+        for name, particle_params in particles.items():
+            E = np.linspace(particle_params['E0_min'], particle_params['E0_max'], 100000) * 1e6 * e  # Umrechnung in Joule
+            dEdx = bethebloch(E, particle_params['m0'], particle_params['z'], params['ne'], params['I'])
+            dEdx_cm = dEdx / (1e6 * e) * 1e-2  # eV/m to MeV/cm
+            linewidth = 1.5 if material != 'BGO' else 3
+            plt.plot(E / (1e6 * e), dEdx_cm, color=colors[name], linestyle=linestyles[material], linewidth=linewidth)
+    
+            if f'{name}' not in legend_particle_labels:
+                legend_particle_handles.append(plt.Line2D([0], [0], color=colors[name], linestyle='-', label=f'{name}'))
+                legend_particle_labels.append(f'{name}')
+    
+        legend_material_handles.append(plt.Line2D([0], [0], color='black', linestyle=linestyles[material], label=material))
+        legend_material_labels.append(material)
+    
+    plt.xscale('log')
+    plt.yscale('log')
+    plt.xlabel('kinetic energy E [MeV]')
+    plt.ylabel('Stopping Power in MeV/cm')
+    plt.title('Energy loss in Si, BGO, Bi, Ge, and O')
+    legend_particle = plt.legend(handles=legend_particle_handles, labels=legend_particle_labels, loc='upper center', title="Particles", frameon=True, fontsize=12*1.5)
+    legend_material = plt.legend(handles=legend_material_handles, labels=legend_material_labels, loc='lower left', title="Materials", frameon=True, fontsize=12*1.5)
+    plt.gca().add_artist(legend_particle)
+    plt.gca().add_artist(legend_particle)
+    
+    plt.grid(True, which="both", ls="--", lw=0.5)
+    
+    ####################################
+    
+    # dE/dx vs. beta * gamma für alle Materialien
+    plt.figure(figsize=(10,6))
+    beta_gamma_range = np.linspace(0.1, 1000, 100000)
+    
+    for material, params in materials.items():
+        for name, particle_params in particles.items():
+            beta_gamma = beta_gamma_range
+            beta = np.sqrt(beta_gamma**2 / (1 + beta_gamma**2))  # Umkehrung von beta*gamma => beta
+            gam = np.sqrt(1 / (1 - beta**2))  # Umkehrung von beta => gamma
+            E = particle_params['m0'] * c**2 * (gam - 1)  # Berechnung der kinetischen Energie aus gamma
+            dEdx = bethebloch(E, particle_params['m0'], particle_params['z'], params['ne'], params['I'])  # Energieverlust pro Distanz
+            dEdx_cm = dEdx / (1e6 * e) * 1e-2  # Umrechnung von eV/m in MeV/cm
+            linewidth = 1.5 if material != 'BGO' else 3
+            plt.plot(beta_gamma, dEdx_cm, color=colors[name], linestyle=linestyles[material], linewidth=linewidth)
+    
+    # Plot-Einstellungen für den zweiten Plot
+    plt.xscale('log')
+    plt.yscale('log')
+    plt.xlabel(r'$\beta \cdot \gamma$')
+    plt.ylabel('Stopping Power in MeV/cm')
+    plt.title('Energy loss in Si, BGO, Bi, Ge, and O (βγ)')
+    legend_particle = plt.legend(handles=legend_particle_handles, labels=legend_particle_labels, loc='upper center', title="Particles", frameon=True, fontsize=12*1.5)
+    legend_material = plt.legend(handles=legend_material_handles, labels=legend_material_labels, loc='upper right', title="Materials", frameon=True, fontsize=12*1.5)
+    plt.gca().add_artist(legend_particle)
+    plt.gca().add_artist(legend_material)
+    plt.grid(True, which="both", ls="--", lw=0.5)
+    
+    plt.show()
+    
+    
+def plot_only_Si():
+    # Plot für alle 3 Teilchen in Silizium
+    plt.figure(figsize=(10,6))
+    colors = {
+        'Proton': 'red',
+        'Helium': 'green',
+        'Muon': (0.5, 0.3, 0.8)  # RGB-Farbdefinition für Muon
+    }
+    # dE/dx vs. Energie (in MeV)
+    for name, params in particles.items():
+        E = np.linspace(params['E0_min'], params['E0_max'], 100000) * 1e6 * e  # Umrechnung in Joule
+        dEdx = bethebloch(E, params['m0'], params['z'], ne_si, I_si)
         dEdx_cm = dEdx / (1e6 * e) * 1e-2  # Umrechnung von eV/m in MeV/cm
-        # Plot für Teilchen mit entsprechender Farbe
+        plt.plot(E / (1e6 * e), dEdx_cm, label=f'{name}', color=colors[name])  # Achsen in MeV
-        linewidth = 2 if material != 'BGO' else 4  # Verdopple die Linienstärke für BGO
+    
-        plt.plot(E / (1e6 * e), dEdx_cm, color=colors[name], linestyle=linestyles[material], linewidth=linewidth)
+    plt.xscale('log')
-
+    plt.yscale('log')
-        # Einträge für die erste Legende (Teilchen), nur für Proton, Helium, und Muon
+    plt.xlabel('kinetic energy E [MeV]')
-        if f'{name}' not in legend_particle_labels:  # Stelle sicher, dass nur einmal hinzugefügt wird
+    plt.ylabel('Stopping Power in MeV/cm')
-            legend_particle_handles.append(plt.Line2D([0], [0], color=colors[name], linestyle='-', label=f'{name}'))
+    plt.title('Energy loss in silicon')
-            legend_particle_labels.append(f'{name}')
+    plt.legend()
-
+    plt.grid(True, which="both", ls="--", lw=0.5)
-    # Material mit schwarzem Linienstil in der zweiten Legende
+    
-    legend_material_handles.append(plt.Line2D([0], [0], color='black', linestyle=linestyles[material], label=material))
+    # dE/dx vs. beta * gamma
-    legend_material_labels.append(material)
+    plt.figure(figsize=(10,6))
-
+    
-# Plot-Einstellungen
+    # Der Bereich für beta * gamma (von 0.1 bis 1000)
-plt.xscale('log')
+    beta_gamma_range = np.linspace(0.1, 1000, 100000)
-plt.yscale('log')
+    
-plt.xlabel('kinetic energy E [MeV]')
+    for name, params in particles.items():
-plt.ylabel('Stopping Power in MeV/cm')
+        # Berechnung von beta und gamma für den Bereich von beta_gamma_range
-plt.title('Energy loss in Si, BGO, Bi, Ge, and O')
+        beta_gamma = beta_gamma_range  # Da beta * gamma bereits gegeben ist
-
-# Erste Legende für die Teilchen (mit Farbkennzeichnung)
-legend_particle = plt.legend(handles=legend_particle_handles, labels=legend_particle_labels, loc='upper center', title="Particles", frameon=True, fontsize=12*1.5)
-plt.gca().add_artist(legend_particle)  # Stelle sicher, dass die erste Legende hinzugefügt wird
-
-# Zweite Legende für Materialien (mit schwarzen Linien)
-legend_material = plt.legend(handles=legend_material_handles, labels=legend_material_labels, loc='lower right', title="Materials", frameon=True, bbox_to_anchor=(0.75, 0.5), fontsize=12*1.5)
-
-plt.grid(True, which="both", ls="--", lw=0.5)
-
-# dE/dx vs. beta * gamma für alle Materialien
-plt.figure(figsize=(10,6))
-
-beta_gamma_range = np.linspace(0.1, 1000, 100000)
-
-for material, params in materials.items():
-    for name, particle_params in particles.items():
-        beta_gamma = beta_gamma_range
         beta = np.sqrt(beta_gamma**2 / (1 + beta_gamma**2))  # Umkehrung von beta*gamma => beta
         gam = np.sqrt(1 / (1 - beta**2))  # Umkehrung von beta => gamma
-        E = particle_params['m0'] * c**2 * (gam - 1)  # Berechnung der kinetischen Energie aus gamma
+        E = params['m0'] * c**2 * (gam - 1)  # Berechnung der kinetischen Energie aus gamma
-        dEdx = -bethebloch(E, particle_params['m0'], particle_params['z'], params['ne'], params['I'])  # Energieverlust pro Distanz
+        dEdx = bethebloch(E, params['m0'], params['z'], ne_si, I_si)  # Energieverlust pro Distanz
         dEdx_cm = dEdx / (1e6 * e) * 1e-2  # Umrechnung von eV/m in MeV/cm
-        # Plot für Teilchen mit entsprechender Farbe
+        plt.plot(beta_gamma, dEdx_cm, label=f'{name}', color=colors[name])  # Achsen in MeV
-        linewidth = 2 if material != 'BGO' else 4  # Verdopple die Linienstärke für BGO
+    
-        plt.plot(beta_gamma, dEdx_cm, color=colors[name], linestyle=linestyles[material], linewidth=linewidth)
+    # Plot-Einstellungen für den zweiten Plot
+    plt.xscale('log')
+    plt.yscale('log')
+    plt.xlabel(r'$\beta \cdot \gamma$')
+    plt.ylabel('Stopping Power in MeV/cm')
+    plt.title('Energy loss in silicon (βγ)')
+    plt.legend()
+    plt.grid(True, which="both", ls="--", lw=0.5)
+    
+    plt.show()
+    
 
-# Plot-Einstellungen für den zweiten Plot
+#CHAOS: A = Si, C= Si, D = BGO, E = Si
-plt.xscale('log')
+
+def calculate_CHAOS(params):
+    E0s = np.linspace(params['E0_min'], params['E0_max'], 1000)
+    dx = 1e-7
+    results = []
+    Si_data = materials['Si']
+    BGO_data = materials['BGO']
+    for E0 in E0s:
+        EA_loss = Eloss(E0, params['m0'], params['z'],  Si_data['ne'], Si_data['I'], dx, Si_data['d'])
+        E = E0 - EA_loss
+        EC_loss = Eloss(E, params['m0'], params['z'],  Si_data['ne'], Si_data['I'], dx, Si_data['d'])
+        E = E - EC_loss
+        ED_loss = Eloss(E, params['m0'], params['z'],  BGO_data['ne'], BGO_data['I'], dx, BGO_data['d'])
+        E = E - ED_loss
+        EE_loss = Eloss(E, params['m0'], params['z'],  Si_data['ne'], Si_data['I'], dx, Si_data['d'])
+        E = E - EE_loss
+        results.append([E0, EA_loss, EC_loss, ED_loss, EE_loss, E])
+    return results
+
+
+results_p = calculate_CHAOS(particles['Proton'])
+results_h = calculate_CHAOS(particles['Helium'])
+results_mu = calculate_CHAOS(particles['Muon'])
+
+def format_frame(results, particle_name):
+    df = pd.DataFrame(results, columns=["E0", "EA", "EC", "ED", "EE", "E_rest"])
+    df = df.round(3)
+    df.to_csv(f'hists_sim/BB-{particle_name}.hist', sep=' ')
+    return df
+
+frame_p = format_frame(results_p, 'Proton')
+frame_h = format_frame(results_h, 'Helium')
+frame_mu = format_frame(results_mu, 'Muon')
+
+print(frame_p)
+print(frame_h)
+print(frame_mu)
+
+plt.plot(frame_p["ED"], frame_p["EC"], 'x', label="proton")
+plt.plot(frame_h["ED"], frame_h["EC"], 'x', label="helium")
+plt.plot(frame_mu["ED"], frame_h["EC"], 'x', label="muon")
+
+plt.xlabel('Energieverlust im BGO (MeV)')
+#plt.ylabel('Verhältnis der Energieverluste (E_si2 / E_si1)')
+plt.ylabel('Energieverluste in (A)')
+plt.title('Energieverluste von Teilchen in Detektoren')
 plt.yscale('log')
-plt.xlabel(r'$\beta \cdot \gamma$')
+plt.ylim(0.01, 100)
-plt.ylabel('Stopping Power in MeV/cm')
+plt.legend()
-plt.title('Energy loss in Si, BGO, Bi, Ge, and O (βγ)')
-
-# Wiederholung der Legenden für beide
-legend_particle = plt.legend(handles=legend_particle_handles, labels=legend_particle_labels, loc='upper right', title="Particles", frameon=True, fontsize=12*1.5)
-plt.gca().add_artist(legend_particle)
-
-legend_material = plt.legend(handles=legend_material_handles, labels=legend_material_labels, loc='lower left', title="Materials", frameon=True, bbox_to_anchor=(0.75, 0.5), fontsize=12*1.5)
-
-plt.grid(True, which="both", ls="--", lw=0.5)
-
 plt.show()
-
+    
-# # Plot für alle 3 Teilchen in Silizium
+    
-# plt.figure(figsize=(10,6))
+    
-
-# colors = {
-#     'Proton': 'red',
-#     'Helium': 'green',
-#     'Muon': (0.5, 0.3, 0.8)  # RGB-Farbdefinition für Muon
-# }
-
-# # dE/dx vs. Energie (in MeV)
-# for name, params in particles.items():
-#     E = np.linspace(params['E0_min'], params['E0_max'], 100000) * 1e6 * e  # Umrechnung in Joule
-#     dEdx = -bethebloch(E, params['m0'], params['z'], ne_si, I_si)
-#     dEdx_cm = dEdx / (1e6 * e) * 1e-2  # Umrechnung von eV/m in MeV/cm
-#     plt.plot(E / (1e6 * e), dEdx_cm, label=f'{name}', color=colors[name])  # Achsen in MeV
-
-# plt.xscale('log')
-# plt.yscale('log')
-# plt.xlabel('kinetic energy E [MeV]')
-# plt.ylabel('Stopping Power in MeV/cm')
-# plt.title('Energy loss in silicon')
-# plt.legend()
-# plt.grid(True, which="both", ls="--", lw=0.5)
-
-# # dE/dx vs. beta * gamma
-# plt.figure(figsize=(10,6))
-
-# # Der Bereich für beta * gamma (von 0.1 bis 1000)
-# beta_gamma_range = np.linspace(0.1, 1000, 100000)
-
-# for name, params in particles.items():
-#     # Berechnung von beta und gamma für den Bereich von beta_gamma_range
-#     beta_gamma = beta_gamma_range  # Da beta * gamma bereits gegeben ist
-#     beta = np.sqrt(beta_gamma**2 / (1 + beta_gamma**2))  # Umkehrung von beta*gamma => beta
-#     gam = np.sqrt(1 / (1 - beta**2))  # Umkehrung von beta => gamma
-#     E = params['m0'] * c**2 * (gam - 1)  # Berechnung der kinetischen Energie aus gamma
-#     dEdx = -bethebloch(E, params['m0'], params['z'], ne_si, I_si)  # Energieverlust pro Distanz
-#     dEdx_cm = dEdx / (1e6 * e) * 1e-2  # Umrechnung von eV/m in MeV/cm
-#     plt.plot(beta_gamma, dEdx_cm, label=f'{name}', color=colors[name])  # Achsen in MeV
-
-# # Plot-Einstellungen für den zweiten Plot
-# plt.xscale('log')
-# plt.yscale('log')
-# plt.xlabel(r'$\beta \cdot \gamma$')
-# plt.ylabel('Stopping Power in MeV/cm')
-# plt.title('Energy loss in silicon (βγ)')
-# plt.legend()
-# plt.grid(True, which="both", ls="--", lw=0.5)
-
-# plt.show()

CHAOSa.py

@@ -38,10 +38,11 @@ parser.add_argument("-cut",type=str, nargs='+',help=("Choose cut definitions for
         "- E123_two: For detector group E123, requires at least two of E1, E2, or E3 to be above threshold."
     )
 )
+parser.add_argument("-fishmask",type=str, nargs='+',help=("Choose mask definitions and apply corresponding conditions:\n"
+        "TBD"))
+
 parser.add_argument("-mask",type=str, nargs='+',help=("Choose mask definitions and apply corresponding conditions:\n"
-        "- HeA: Applies He-mask using A/E123 vs D (fishplot).\n"
+        "- PlowE: Applies mask for lower energetic Protons in min(A,C)-D plot.\n"))
-        "- HeC: Applies He-mask using C/E123 vs D (fishplot).\n"
-        "- He1: Applies He-mask for min(A,C)-D plot"))
 
 parser.add_argument("-deltaT", action='store_true', help = "saves deltatimes between events")
 
@@ -74,12 +75,12 @@ from histframe import get_T
 from histframe import CmV
 from histframe import CmV_v
 from histframe import isCut
-from histframe import isMask
-from histframe import isMask_other
 from histframe import isMask_1D
 from histframe import get_MeV_mean
 from histframe import get_MeV
 
+from masks import *
+
 
 names = hf.names
 ch = hf.ch
@@ -246,8 +247,8 @@ def main():
         generate_2d_hist_logxy()
 
 
-def conditions(l,args): 
+def conditions(l,T,args): 
-    return (not args.cut or isCut(l,args.cut)) and (not args.trig or isTrigger(l)) and (not args.xray or not isBGO(l))
+    return (not args.mask or isMask_other(l,T,args.mask)) and (not args.cut or isCut(l,args.cut)) and (not args.trig or isTrigger(l)) and (not args.xray or not isBGO(l))
 
 # def isTrigger(l):
 #     return all( any(CmV(l, ch[sub_ch]) > thr[sub_ch] for sub_ch in chans) for chans in trigchan)
@@ -338,7 +339,7 @@ def generate_hist():
                 if l[0] == 'H':
                     timestamp = int(l[1])
                     T = get_T(float(l[8]))
-                if l[0] == 'EI' and len(l)>3*18 and conditions(l,args) and isMask_other(l,T,mask) and isFloat(args,timestamp):
+                if l[0] == 'EI' and len(l)>3*18 and conditions(l,T,args) and isMask_other(l,T,mask) and isFloat(args,timestamp):
                     for i,chan in enumerate(chans):
                         b = value_HL(l,chan,u) 
                         #b = value
@@ -375,7 +376,7 @@ def generate_hist_D_MeV():
                 if l[0] == 'H':
                     timestamp = int(l[1])
                     T = get_T(float(l[8]))
-                if l[0] == 'EI' and len(l)>3*18 and conditions(l,args) and isFloat(args,timestamp):
+                if l[0] == 'EI' and len(l)>3*18 and conditions(l,T,args) and isFloat(args,timestamp):
                     #T=15
                     b = value_HL(l,'D1',u)+value_HL(l,'D2',u) 
                     bx = get_MeV(b, T)
@@ -413,7 +414,7 @@ def generate_hist_D_MeV_log():
                 if l[0] == 'H':
                     timestamp = int(l[1])
                     T = get_T(float(l[8]))
-                if l[0] == 'EI' and len(l)>3*18 and conditions(l,args):
+                if l[0] == 'EI' and len(l)>3*18 and conditions(l,T,args):
                     #T=15
                     b = value_HL(l,'D1',u)+value_HL(l,'D2',u) 
                     bx = get_MeV(b, T)
@@ -431,9 +432,9 @@ def generate_hist_D_MeV_log():
 def generate_hist_xlog():
     minX = 1
     maxX = 3500
-    resX = 100*resV
+    maxX = 100000
+    resX = 1000*resV
     bin_x = int((maxX-minX)/resX)
-    
     edges_x = np.logspace(np.log10(minX), np.log10(maxX), bin_x)
     
     hist = np.zeros((bin_x,len(chans)+1))
@@ -452,13 +453,13 @@ def generate_hist_xlog():
                 if l[0] == 'H':
                     timestamp = int(l[1])
                     T = get_T(float(l[8]))
-                if l[0] == 'EI' and len(l)>3*18 and conditions(l,args) and isMask_other(l,T,mask):
+                if l[0] == 'EI' and len(l)>3*18 and conditions(l,T,args) and isMask_other(l,T,mask):
                     for i,chan in enumerate(chans):
                         #b = value_HL(l,chan,u) 
                         x_expr, x_op, x_next_op, x_next_channel = parse_channel_expression(chan)
                         try: b = get_channel_value(l, u,x_expr, x_op, x_next_op, x_next_channel,T)
-                        except: b = 0
+                        except: b = maxX+1
-                        if minV <= b <= maxX:
+                        if minX<= b <= maxX:
                             xb = np.digitize(b, edges_x)-1
                             hist[xb,i+1] += 1                                                           
                                 
@@ -494,7 +495,7 @@ def generate_2d_hist():
             if l[0] == 'H':
                 timestamp = int(l[1])
                 T = get_T(float(l[8]))
-            if l[0] == 'EI' and len(l)>3*18 and conditions(l,args) and isFloat(args,timestamp) and T != None:  
+            if l[0] == 'EI' and len(l)>3*18 and conditions(l,T,args) and isFloat(args,timestamp) and T != None:  
                 if int(int(l[4].split(':')[0])/3) > deltacut:
                     bx = get_channel_value(l, u,x_expr, x_op, x_next_op, x_next_channel,T)
                     by = get_channel_value(l, u, y_expr, y_op, y_next_op, y_next_channel,T)
@@ -548,7 +549,7 @@ def generate_2d_hist_logy():
             if l[0] == 'H':
                 timestamp = int(l[1])
                 T = get_T(float(l[8]))
-            if l[0] == 'EI' and len(l)>3*18 and conditions(l,args) and isFloat(args,timestamp) and T != None:   
+            if l[0] == 'EI' and len(l)>3*18 and conditions(l,T,args) and isFloat(args,timestamp) and T != None:   
                 if int(int(l[4].split(':')[0])/3) > deltacut:
                     bx = get_channel_value(l,u, x_expr, x_op, x_next_op, x_next_channel,T)
                     by = get_channel_value(l,u, y_expr, y_op, y_next_op, y_next_channel,T)
@@ -614,7 +615,7 @@ def generate_2d_hist_logxy():
             if l[0] == 'H':
                 timestamp = int(l[1])
                 T = get_T(float(l[8]))
-            if l[0] == 'EI' and len(l)>3*18 and conditions(l,args) and isFloat(args,timestamp) and T != None:   
+            if l[0] == 'EI' and len(l)>3*18 and conditions(l,T,args) and isFloat(args,timestamp) and T != None:   
                 if int(int(l[4].split(':')[0])/3) > deltacut:
                     bx = get_channel_value(l,u, x_expr, x_op, x_next_op, x_next_channel,T)
                     by = get_channel_value(l,u, y_expr, y_op, y_next_op, y_next_channel,T)

CHAOScustom.py

@@ -15,6 +15,8 @@ from histframe import CmV
 from histframe import isCut
 from histframe import isMask
 
+from masks import isMask_other
+
 
 name = hf.names
 ch = hf.ch
@@ -28,11 +30,12 @@ resV = hf.resV
 
 pol = hf.pol
 
-file = "2024-10-02-flight-2.EI"
+# file = "2024-10-02-flight-2.EI"
-# file = "chaos_sd_float.EI"
+file = "chaos_sd_float.EI"
-# file = "data_sim/proton.EI"
+file = "data_sim/proton.EI"
-# file = "data_sim/e-.EI"
+#file = "data_sim/e-.EI"
-# file = "data_sim/mu.EI"file = "data_sim/he.EI"
+#file = "data_sim/mu.EI"
+#file = "data_sim/he.EI"
 
 
 def do_hist(cuts,trigchans,val):
@@ -60,7 +63,7 @@ def do_hist(cuts,trigchans,val):
             if l[0] == 'H':
                 timestamp = int(l[1])
                 T = get_T(float(l[8]))
-            if l[0] == 'EI' and len(l)>3*18 and isCut(l,cuts) and isTrigger(l,trigchans) and T != None:
+            if l[0] == 'EI' and len(l)>3*18 and isCut(l,cuts) and isMask_other(l,T,masks) and isTrigger(l,trigchans) and T != None:
                 if int(int(l[4].split(':')[0])/3) > deltacut:
                     if val=='fish':
                         bx,by = get_values_fish(l,T)
@@ -68,7 +71,7 @@ def do_hist(cuts,trigchans,val):
                         bx,by = get_values_other(l,T)
                     else:
                         bx,by = get_values(l,T)
-                    if minX <= bx <= maxX and minY <= by <= maxY and isMask(bx,by,masks):
+                    if minX <= bx <= maxX and minY <= by <= maxY:
                         bxx = np.digitize(bx, edges_x) - 1
                         #bxx = int((bx-minX)/resX)
                         byy = np.digitize(by, edges_y)
@@ -84,15 +87,15 @@ def isTrigger(l, trigchans):
     return (all(any(CmV(l, ch[sub_ch]) > thr[sub_ch] for sub_ch in chan) for chan in trigchans))
 
 def get_values_fish(l,T):
-    y = min(value_HL(l,"A",u),value_HL(l,"C",u))/max(value(l,"E1",u),value(l,"E2",u),value(l,"E3",u)) 
+    #y = min(value_HL(l,"A",u),value_HL(l,"C",u))/max(value(l,"E1",u),value(l,"E2",u),value(l,"E3",u)) 
-    #y = value_HL(l,"A",u)/max(value(l,"E1",u),value(l,"E2",u),value(l,"E3",u)) 
+    y = value_HL(l,"C",u)/max(value(l,"E1",u),value(l,"E2",u),value(l,"E3",u)) 
     x = value_HL(l,"D1",u)+value_HL(l,"D2",u)
     x2 = get_MeV(x,T)
     return x2,y
 
 def get_values_other(l,T):
-    y = min(value_HL(l,"A",u),value_HL(l,"C",u))/1000
+    #y = min(value_HL(l,"A",u),value_HL(l,"C",u))/1000
-    #y = value_HL(l,"A",u)/1000
+    y = value_HL(l,"C",u)/1000
     x = value_HL(l,"D1",u)+value_HL(l,"D2",u)
     x2 = get_MeV(x,T)
     return x2,y
@@ -110,23 +113,24 @@ def get_values(l,T):
 
 deltacut = 0
 
-allcuts = [['A','B','C','AeqC','D'],['A','nB','C','AeqC','D']]
+allcuts = [['A','C','AeqC','D'],['A','B','C','AeqC','D'],['A','nB','C','AeqC','D'],['A','rB','C','AeqC','D']]
-allcuts = [['A','B','C','D'],['A','nB','C','D']]
+#allcuts = [['A','Bsim','C','AeqC','D'],['A','nBsim','C','AeqC','D']]
 #allcuts = [['A','D'],['A','B','D'],['A','nB','D']]
-#allcuts = [['A','C','AeqC','D']]
+allcuts = [['A','nB','C','AeqC','D']]
 
-trig = ['A1H','C1H']
+#trig = ['A1H','C1H']
 trig=[]
 trigchans = [[i for sub_t in t.split("v") for i, n in enumerate(hf.names) if sub_t == n] for t in trig]
 
-#masks = ['He1']
 masks=[]
+#masks = ['He1']
+#masks=["PlowE"]
 
 allfish = True
 allother = True
 
 folder = "hists_fish"
-#folder = "hists_sim"
+folder = "hists_sim"
 
 minX = 2
 maxX = 10000
@@ -138,25 +142,23 @@ resX = resV/0.00362/polynom(15,*pol)
 if allother:
     for cut in allcuts:
         filename = "chaosB43_Cut"+ "".join(cut)+"~min(A,C)-D.2dhist2log"
+        filename = "chaosB43_Cut"+ "".join(cut)+"~C-D.2dhist2log"
+        #filename = "chaosB43_Cut"+ "".join(cut)+'_mask'+"".join(masks)+"~min(A,C)-D.2dhist2log"
         #filename = "chaos_Cut"+ "".join(cut)+"_trig-"+ "".join(trig)+"~min(A,C)-D.2dhist2log"
-        #filename = "chaosB43_Cut"+ "".join(cut)+"~A-D.2dhist2log"
+        filename = "simProton_Cut"+ "".join(cut)+"~C-D.2dhist2log"
-        #filename = "simHe-chaos_Cut"+ "".join(cut)+"~min(A,C)-D.2dhist2log"
         do_hist(cut,trigchans,'other') 
         
 if allfish:
     for cut in allcuts:
         cut.append('E123')
         filename = "chaosB43_Cut"+ "".join(cut)+"~min(A,C)dE-D.2dhist2log"
+        filename = "chaosB43_Cut"+ "".join(cut)+"~CdE-D.2dhist2log"
+        #filename = "chaosB43_Cut"+ "".join(cut)+'_mask'+"".join(masks)+"~min(A,C)dE-D.2dhist2log"
         #filename = "chaos_Cut"+ "".join(cut)+"_trig-"+ "".join(trig)+"~min(A,C)dE-D.2dhist2log"
-        #filename = "chaosB43_Cut"+ "".join(cut)+"~AdE-D.2dhist2log"
+        filename = "simProton_Cut"+ "".join(cut)+"~CdE-D.2dhist2log"
-        #filename = "simHe-chaos_Cut"+ "".join(cut)+"~min(A,C)dE-D.2dhist2log"
         do_hist(cut,trigchans,'fish')
 
 
-
-        
-        
-
 #cuts = ['A','hB','C','AeqC','D','E123']
 
 #filename = "chaos35_Cut"+ "".join(cuts)+"~min(A,C)dE-D.2dhist2log"

CHAOSx.py

@@ -278,9 +278,9 @@ def generate_deltaTimes():
 ################ 1D HISTOGRAM ############################################################
 
 def generate_hist():
-    minV = 0
+    minV = -100
-    maxV = 30
+    maxV = 10000
-    resV = 0.1
+    #resV = 10*resV
     bin = int((maxV-minV)/resV)
     
     hist = np.zeros((bin+1,18+1))

histframe.py

@@ -156,7 +156,7 @@ def isCut(l, cut):
         
     if "B" in cut or "hB" in cut:
         b = value(l,"B",u)
-        if b < 17: return False
+        if b < 43: return False
         #if b > 106.38: return False
         
     elif "nB" in cut or "lB" in cut:
@@ -167,6 +167,15 @@ def isCut(l, cut):
         b = value(l,"B",u)
         if b < 17 or b > 43: return False
         
+    if "Bsim" in cut:
+        b = value(l,"B",u)
+        if b < 227*0.65: return False
+        #if b > 106.38: return False
+        
+    elif "nBsim" in cut:
+        b = value(l,"B",u)
+        if b > 227*0.65: return False
+        
     if "A" in cut:
         a1 = value_HL(l,"A1",u)
         a2 = value_HL(l,"A2",u)

masks.py

@@ -0,0 +1,108 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Feb 19 10:53:26 2025
+
+@author: ava
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+
+def linfkt(p1,p2):
+    return lambda x: (p2[1]-p1[1])/(p2[0]-p1[0]) * (x-p1[0]) +p1[1]
+
+def logloglinfkt(p1,p2):
+    m=np.log(p2[1]/p1[1])/np.log(p2[0]/p1[0])
+    b=p1[1]/p1[0]**m
+    return lambda x: b*x**m
+
+def mask(box): #from top left to bottom left clockwise, only shapes with 4 sides for now
+    if len(box)==4:
+        p1=box[0]
+        p2=box[1]
+        p3=box[2]
+        p4=box[3]
+        g=logloglinfkt(p1,p2)
+        f=logloglinfkt(p3,p4)
+        j=logloglinfkt(p1[::-1], p4[::-1])
+        i=logloglinfkt(p2[::-1], p3[::-1])
+        return (g,f,j,i)   
+    
+def drawbox(ax,box,color="b"):
+    # dbox=box(lambda x:x,lambda y:y,p1,p2,p3,p4)
+    if len(box)==4:      
+        p1=box[0]
+        p2=box[1]
+        p3=box[2]
+        p4=box[3]
+        gx=np.linspace(p1[0],p2[0],1000)
+        fx=np.linspace(p3[0],p4[0],1000)
+        iy=np.linspace(p2[1],p3[1],1000)
+        jy=np.linspace(p4[1],p1[1],1000)
+        ax.plot(gx,logloglinfkt(p1,p2)(gx),color='white', linewidth=3)
+        ax.plot(fx,logloglinfkt(p3,p4)(fx),color='white', linewidth=3)
+        ax.plot(logloglinfkt(p1[::-1],p4[::-1])(jy),jy,color='white', linewidth=3)
+        ax.plot(logloglinfkt(p2[::-1],p3[::-1])(iy),iy,color='white', linewidth=3)
+        ax.plot(gx,logloglinfkt(p1,p2)(gx),color=color, linewidth=1.5)
+        ax.plot(fx,logloglinfkt(p3,p4)(fx),color=color, linewidth=1.5)
+        ax.plot(logloglinfkt(p1[::-1],p4[::-1])(jy),jy,color=color, linewidth=1.5)
+        ax.plot(logloglinfkt(p2[::-1],p3[::-1])(iy),iy,color=color, linewidth=1.5)
+        
+#niederenergetische Protonen
+PlowE1=(9,1.1)
+PlowE2=(67,0.7)
+PlowE3=(39,0.4)
+PlowE4=(7.5,0.8)
+box_PlowE=(PlowE1,PlowE2,PlowE3,PlowE4)
+
+# print(PlowE2[::-1])
+# print(PlowE3[::-1])
+#print(mask(box_PlowE))
+
+#MIP
+mip1=(12.5,0.14)
+mip2=(28,0.14)
+mip3=(28,0.05)
+mip4=(12.5,0.05)
+box_mip=(mip1,mip2,mip3,mip4)
+
+
+from histframe import get_values_other       
+
+def isMask_other(l,T,masks):
+    x,y = get_values_other(l, T)
+    if "PlowE" in masks:
+        g,f,j,i = mask(box_PlowE)
+        #if f(x)<y<g(x) and i(y) > x > j(y): return True
+        if y<f(x) or y > g(x) or x > i(y) or x < j(y) : return False
+        
+    if "MIP" in masks:
+        g,f,j,i = mask(box_mip)
+        if y<f(x) or y > g(x) or x > i(y) or x < j(y) : return False
+        
+    return True
+
+from histframe import get_values_fish           
+def isMask_fish(l, T,masks):
+    if "PlowE" in masks:
+        x,y = get_values_fish(l, T)
+    return True
+
+######################################################################
+def count_points_in_box(box, x_edges, y_edges, c):
+    p1, p2, p3, p4 = box
+    g,f,j,i = mask(box)
+
+    x_centers = (np.array(x_edges[:-1]) + np.array(x_edges[1:])) / 2
+    y_centers = (np.array(y_edges[:-1]) + np.array(y_edges[1:])) / 2
+    total_count = 0
+    for ix, x_center in enumerate(x_centers):
+        for iy, y_center in enumerate(y_centers):
+            if (p1[0] <= x_center <= p2[0] and f(x_center) <= y_center <= g(x_center)):
+                total_count += c[iy, ix]
+
+    return total_count
+
+                

plot3d.py

@@ -0,0 +1,245 @@
+'''
+Autor: Ava Pohley
+'''
+
+import argparse
+import sys
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+import pandas as pd
+import numpy as np
+import os.path
+from matplotlib.colors import ListedColormap, Normalize
+from scipy.optimize import curve_fit
+from matplotlib.ticker import FuncFormatter
+from matplotlib.gridspec import GridSpec
+
+from masks import *
+ 
+parser = argparse.ArgumentParser(description="This script plots histogramms of given data ending on '.2dhist'")
+parser.add_argument("-v", "--verbose", choices = ['0','1','2'], help = "when zero only emit serious error messages; when one also emit warnings, when two emit informational messages")
+parser.add_argument("filepath", type = str, help="path/file were data is stored in")
+parser.add_argument("-t", "--title", type = str, help="title to use in plot; if not given a standard title will be chosen")
+parser.add_argument("-e", "--export", type = str, nargs='+', help = "first: filepath were figure shall be exported to (if f, filepath of histdata will be chosen with filepath_plot) \n second:format (png,pdf,...)")
+parser.add_argument("-l", "--limits", type=float, nargs='+',help = "limits of plot in order xmin,xmax, ymin,ymax,")
+parser.add_argument("-d", "--diagonal",action='store_true', help = "draws black line x=y")
+parser.add_argument("-c", "--channels", type=str, nargs='+',help = "only use for .hist files (1d). letters of chanels to be plot; if not given, all channels will be plotted")
+parser.add_argument("-p", "--photoelectrons", action='store_true', help = "converts data from keV{Si} to amount of photoelectrons")
+parser.add_argument("-k", "--kontour", action='store_true', help = "plots kontour instead of 2d hist")
+parser.add_argument("-f", "--files", type=str, nargs='+',help = "only use for .hist files (1d). add files from wich data shall be also plotted")
+parser.add_argument("-n", "--names", type=str, nargs='+',help = "only use for -f. Specifies labels of channes in order of appearence")
+
+
+args = parser.parse_args()
+
+if args.verbose == None:
+    verbosity = 1
+else:
+    verbosity = args.verbose
+    
+if verbosity == 2:
+    sys.stderr.write("--verbostiy = '" + str(args.verbose)+"'\n")
+    sys.stderr.write("--filepath = '" + str(args.filepath)+"'\n")
+    sys.stderr.write("--title = '" + str(args.title)+"'\n")
+    sys.stderr.write("--export = '" + str(args.export) +"'\n")
+    sys.stderr.write("--limits = '" + str(args.limits) +"'\n")
+    sys.stderr.write("--channels = '" + str(args.channels)+"'\n")
+    sys.stderr.write("--files = '" + str(args.files)+"'\n")
+
+   
+file, end = args.filepath.split('.')
+
+if args.photoelectrons:
+    f = 1/0.00362
+else: 
+    f = 1.0
+
+def main():
+    if end.startswith('2'):
+        do_2dhist()
+    else: print("skript not usable for filetyp"+end)
+    
+###############################################################################
+        
+def do_2dhist():
+    fs = 25  # Schriftgröße für Labels und Achsen
+    try:
+        y, x = file.split('~')[-1].split('-')
+    except ValueError:
+        xlab, ylab = 'unknown', 'unknown'
+    
+    tit = args.title if args.title else f'Histplot2D, "{file.split("/")[-1].split("~")[0]}"'
+    
+    # Erstelle eine Figur mit GridSpec
+    fig = plt.figure(figsize=(15, 10))
+    spec = GridSpec(nrows=2, ncols=2, width_ratios=[4, 0.1], height_ratios=[6, 1], figure=fig)  # Layout
+    spec = GridSpec(nrows=1, ncols=2, width_ratios=[4, 0.1], figure=fig)  # Hauptplot, Colorbar, Summenplot
+
+   
+    
+    # Hauptplot (oben links)
+    ax_main = fig.add_subplot(spec[0, 0])
+    ax_main.set_xscale('log')
+    
+    # Colorbar (oben rechts)
+    #cbar_ax = fig.add_subplot(spec[0, 1])
+    
+    # Summenplot (unten links)
+    #ax_sum = fig.add_subplot(spec[1, 0], sharex=ax_main)
+    #ax_sum = fig.add_subplot(spec[0, 1], sharey=ax_main)
+    cbar_ax = fig.add_subplot(spec[0, 1])
+
+    
+    # Daten lesen
+    data = read_data(file, end)
+    energy = data.columns[0]
+    x_edges = data[energy].tolist()
+    y_edges = data.columns[1:].astype(float).tolist()
+    cmap = mpl.colormaps['inferno']
+    cmap = mpl.colormaps['viridis']
+    #cmap = mpl.colormaps['magma']
+    cmap.set_under('white')
+    c = np.array(data)[:, 1:].T
+
+    # Hauptplot: 2D-Histogramm
+    #cb = ax_main.pcolormesh(x_edges, y_edges, c, cmap=cmap, shading='nearest', vmin=1,vmax=100, norm='log')
+    cb = ax_main.pcolormesh(x_edges, y_edges, c, cmap=cmap, shading='nearest', vmin=1, norm='log')
+    
+    ax_main.set_yscale('log')  # Logarithmische Skala für y
+    
+    
+    
+###############   DRAW  BOXES    ##############################################################################
+    #drawbox(ax_main,box_PlowE)    
+    #drawbox(ax_main,box_mip)  
+    
+    #counts_in_box = count_points_in_box(box_mip, x_edges, y_edges, c)
+
+    #ax_main.plot(0, 0, 'b-', label=f'MIP box (Counts: {int(counts_in_box)})')
+    #ax_main.legend(fontsize=fs)
+    
+    
+    # Summenberechnung für den unteren Plot
+    x_sum = c[:, :].sum(axis=0)
+    y_sum = c[:, :].sum(axis=1)
+    
+    if 'd' in y:
+        ylab = y.split('d')[0] + ' / ' + y.split('d')[1]
+    else: ylab = f"{y} in {energy}"
+    xlab = f"{x} in {energy}"
+    
+    # Labels und Achsen für Hauptplot
+    ax_main.set_xlabel("", fontsize=fs)  # Kein Label unten, nur oben
+    ax_main.set_ylabel(ylab, fontsize=fs + 5)
+    #ax_main.tick_params(axis='x', labelsize=fs)
+    #ax_main.tick_params(axis='y', labelsize=fs)
+    ax_main.tick_params(axis='x', labelsize=fs,size=7, width=1.5)
+    ax_main.tick_params(axis='y', labelsize=fs,size=7, width=1.5)
+    ax_main.set_title(tit, fontsize=fs + 5)
+    
+    # Colorbar
+    cbar = fig.colorbar(cb, cax=cbar_ax, orientation='vertical')
+    cbar.set_label('counts', fontsize=fs + 5)
+    #cbar.ax.tick_params(labelsize=fs)
+    cbar.ax.tick_params(labelsize=fs,size=7, width=1.5)
+    ax_main.set_xlabel(xlab, fontsize=fs + 5)
+    
+    #Summenplot (unterhalb des Hauptplots)
+    # ax_sum.step(x_edges[:len(x_sum)], x_sum, color='blue')
+    # ax_sum.set_xlabel(xlab, fontsize=fs + 5)  # Gemeinsames x-Label
+    # ax_sum.set_ylabel("counts", fontsize=fs + 5)
+    # ax_sum.tick_params(axis='x', labelsize=fs)
+    # ax_sum.tick_params(axis='y', labelsize=fs)
+    # #ax_sum.set_yscale('log')  # Logarithmische Skala für die Summen
+    
+    
+    # ax_sum.step(y_sum, y_edges[:len(y_sum)], color='blue')  # Zeilensummen gegen y-Werte
+    # ax_sum.set_xlabel("counts", fontsize=fs + 5)
+    # #ax_sum.tick_params(axis='x', labelsize=fs)
+    # #ax_sum.tick_params(axis='y', labelsize=fs)
+    # ax_sum.set_xscale('log')  # Logarithmische Skala für die Summen
+    # ax_sum.tick_params(axis='y', left=False, labelleft=False)
+    # ax_sum.set_xlim(1, 600)
+    
+    
+    AeqC = False
+    if AeqC:
+        a=2
+        ax_main.axhline(y=a, color='r', linestyle = '-', lw=1.5, label='A = '+str(a)+'C  |  C = '+str(a)+'A',)
+        ax_main.axhline(y=1/a, color='r', linestyle = '-', lw=1.5,)
+        a=2
+        ax_sum.axhline(y=a, color='r', linestyle = '-', lw=1.5,)
+        ax_sum.axhline(y=1/a, color='r', linestyle = '-', lw=1.5)
+        ax_main.legend(fontsize=fs, loc="upper left",ncol=2)
+    
+    # Achsenlimits setzen
+    limits = args.limits if args.limits else [x_edges[0], x_edges[-1], y_edges[0], y_edges[-1]]
+    ax_main.set_xlim(limits[0], limits[1])
+    ax_main.set_ylim(limits[2], limits[3])
+    #ax_sum.set_xlim(limits[0], limits[1])  # Gleiche x-Achse wie der Hauptplot
+    
+    # Layout optimieren
+    plt.tight_layout(rect=[0, 0, 0.95, 1])  # Platz für Colorbar rechts
+    save_or_show_plot([])
+
+
+###############################################################################    
+    
+def read_data(file, end):
+    return pd.read_csv(f"{file}.{end}", sep='\s+', skiprows=0)
+ 
+    
+###############################################################################
+
+def read_Itime(file, end):    
+    f = file + '.Itime'
+    Itime = None
+    while len(f) > 0:
+        if os.path.isfile(f+'.Itime'):
+            Itime = int(pd.read_csv(f+'.Itime', sep='\s+', skiprows=0, names = ['time'])['time'][0])
+            f=''
+        else:
+            f=f[:-1]
+    return Itime
+
+    
+###############################################################################
+    
+def save_or_show_plot(channels):
+    if args.export:
+        if args.export[0] != 'f':
+            name = args.export[0] 
+        else:
+            name = f'plots/{end}_{file.split("/")[-1]}'
+            if channels:
+                name += '_ch-' + ''.join(channels)
+        plt.tight_layout()
+        plt.savefig(f"{name}.{args.export[1] if len(args.export) > 1 else 'png'}")
+        print(f'saved as {name}')
+    plt.show()  
+   
+    
+################################
+def cutA1(A1):
+    arr=[]
+    for val in A1:
+        if val<8: arr.append(val+11)
+        else: arr.append((A1-8)/0.008)
+    return arr
+    
+def cutA2(A1):
+    arr=[]
+    for val in A1:
+        if val<8: arr.append((A1-7)/0.15)
+        else: arr.append(0.008*A1+11)
+    return arr
+###############################################################################
+
+def constant(x, c):
+    #return np.full_like(x,c)
+    return c
+
+###############################################################################
+
+if __name__ == '__main__':
+    main()

plot_simulations.py

@@ -0,0 +1,233 @@
+'''
+Autor: Ava Pohley
+'''
+
+import argparse
+import sys
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+import pandas as pd
+import numpy as np
+import os.path
+from matplotlib.colors import ListedColormap, Normalize
+from scipy.optimize import curve_fit
+from matplotlib.ticker import FuncFormatter
+from matplotlib.gridspec import GridSpec
+from matplotlib.colors import LogNorm
+ 
+parser = argparse.ArgumentParser(description="This script plots histogramms of given data ending on '.2dhist'")
+parser.add_argument("-v", "--verbose", choices = ['0','1','2'], help = "when zero only emit serious error messages; when one also emit warnings, when two emit informational messages")
+parser.add_argument("filepath", type = str, help="path/file were data is stored in")
+parser.add_argument("-t", "--title", type = str, help="title to use in plot; if not given a standard title will be chosen")
+parser.add_argument("-e", "--export", type = str, nargs='+', help = "first: filepath were figure shall be exported to (if f, filepath of histdata will be chosen with filepath_plot) \n second:format (png,pdf,...)")
+parser.add_argument("-l", "--limits", type=float, nargs='+',help = "limits of plot in order xmin,xmax, ymin,ymax,")
+parser.add_argument("-d", "--diagonal",action='store_true', help = "draws black line x=y")
+parser.add_argument("-c", "--channels", type=str, nargs='+',help = "only use for .hist files (1d). letters of chanels to be plot; if not given, all channels will be plotted")
+parser.add_argument("-p", "--photoelectrons", action='store_true', help = "converts data from keV{Si} to amount of photoelectrons")
+parser.add_argument("-k", "--kontour", action='store_true', help = "plots kontour instead of 2d hist")
+parser.add_argument("-f", "--files", type=str, nargs='+',help = "only use for .hist files (1d). add files from wich data shall be also plotted")
+parser.add_argument("-n", "--names", type=str, nargs='+',help = "only use for -f. Specifies labels of channes in order of appearence")
+
+
+args = parser.parse_args()
+
+if args.verbose == None:
+    verbosity = 1
+else:
+    verbosity = args.verbose
+    
+if verbosity == 2:
+    sys.stderr.write("--verbostiy = '" + str(args.verbose)+"'\n")
+    sys.stderr.write("--filepath = '" + str(args.filepath)+"'\n")
+    sys.stderr.write("--title = '" + str(args.title)+"'\n")
+    sys.stderr.write("--export = '" + str(args.export) +"'\n")
+    sys.stderr.write("--limits = '" + str(args.limits) +"'\n")
+    sys.stderr.write("--channels = '" + str(args.channels)+"'\n")
+    sys.stderr.write("--files = '" + str(args.files)+"'\n")
+
+   
+file, end = args.filepath.split('.')
+
+if args.photoelectrons:
+    f = 1/0.00362
+else: 
+    f = 1.0
+
+def main():
+    if end.startswith('2'):
+        do_2dhist()
+    else: print("skript not usable for filetyp"+end)
+    
+###############################################################################
+        
+def do_2dhist():
+    fs = 25  # Schriftgröße für Labels und Achsen
+    try:
+        y, x = file.split('~')[-1].split('-')
+    except ValueError:
+        xlab, ylab = 'unknown', 'unknown'
+    
+    tit = args.title if args.title else f'Histplot2D, "{file.split("/")[-1].split("~")[0]}"'
+    
+    # Erstelle eine Figur mit GridSpec
+    fig = plt.figure(figsize=(15, 10))
+    spec = GridSpec(nrows=2, ncols=2, width_ratios=[4, 0.1], height_ratios=[6, 1], figure=fig)  # Layout
+    #spec = GridSpec(nrows=1, ncols=3, width_ratios=[4, 0.8, 0.1], figure=fig)  # Hauptplot, Colorbar, Summenplot
+    spec = GridSpec(nrows=1, ncols=2, width_ratios=[4, 0.1], figure=fig)  # Layout
+   
+    
+    # Hauptplot (oben links)
+    ax_main = fig.add_subplot(spec[0, 0])
+    ax_main.set_xscale('log')
+    
+    # Colorbar (oben rechts)
+    #cbar_ax = fig.add_subplot(spec[0, 1])
+    
+    # Summenplot (unten links)
+    #ax_sum = fig.add_subplot(spec[1, 0], sharex=ax_main)
+    #ax_sum = fig.add_subplot(spec[0, 1], sharey=ax_main)
+    cbar_ax = fig.add_subplot(spec[0, 1])
+
+    
+    # Daten lesen
+    data = read_data(file, end)
+    energy = data.columns[0]
+    x_edges = data[energy].tolist()
+    y_edges = data.columns[1:].astype(float).tolist()
+    cmap = mpl.colormaps['inferno']
+    cmap = mpl.colormaps['viridis']
+    cmap.set_under('white')
+    c = np.array(data)[:, 1:].T
+
+    # Hauptplot: 2D-Histogramm
+    #cb = ax_main.pcolormesh(x_edges, y_edges, c, cmap=cmap, shading='nearest', vmin=1,vmax=100, norm='log')
+    cb = ax_main.pcolormesh(x_edges, y_edges, c, cmap=cmap, shading='nearest', vmin=1, norm='log')
+    
+    ax_main.set_yscale('log')  # Logarithmische Skala für y
+    
+    #frame_p = pd.read_csv('hists_sim/BB-Proton.hist', sep=' ')
+    #sc_p = ax_main.scatter(frame_p["ED"][1:], frame_p["EC"][1:], c=frame_p["E0"][1:], cmap='inferno', norm=LogNorm(), marker='o')
+    
+    # Summenberechnung für den unteren Plot
+    x_sum = c[:, :].sum(axis=0)
+    y_sum = c[:, :].sum(axis=1)
+    
+    if 'd' in y:
+        ylab = y.split('d')[0] + ' / ' + y.split('d')[1]
+    else: ylab = f"{y} in {energy}"
+    xlab = f"{x} in {energy}"
+    
+    # Labels und Achsen für Hauptplot
+    ax_main.set_xlabel("", fontsize=fs)  # Kein Label unten, nur oben
+    ax_main.set_ylabel(ylab, fontsize=fs + 5)
+    #ax_main.tick_params(axis='x', labelsize=fs)
+    #ax_main.tick_params(axis='y', labelsize=fs)
+    ax_main.tick_params(axis='x', labelsize=fs,size=7, width=1.5)
+    ax_main.tick_params(axis='y', labelsize=fs,size=7, width=1.5)
+    ax_main.set_title(tit, fontsize=fs + 5)
+    
+    # Colorbar
+    cbar = fig.colorbar(cb, cax=cbar_ax, orientation='vertical')
+    cbar.set_label('counts', fontsize=fs + 5)
+    #cbar.ax.tick_params(labelsize=fs)
+    cbar.ax.tick_params(labelsize=fs,size=7, width=1.5)
+    
+    # Summenplot (unterhalb des Hauptplots)
+    # ax_sum.step(x_edges[:len(x_sum)], x_sum, color='blue')
+    # ax_sum.set_xlabel(xlab, fontsize=fs + 5)  # Gemeinsames x-Label
+    # ax_sum.set_ylabel("counts", fontsize=fs + 5)
+    # ax_sum.tick_params(axis='x', labelsize=fs)
+    # ax_sum.tick_params(axis='y', labelsize=fs)
+    # #ax_sum.set_yscale('log')  # Logarithmische Skala für die Summen
+    
+    ax_main.set_xlabel(xlab, fontsize=fs + 5)
+    # ax_sum.step(y_sum, y_edges[:len(y_sum)], color='blue')  # Zeilensummen gegen y-Werte
+    # ax_sum.set_xlabel("counts", fontsize=fs + 5)
+    # #ax_sum.tick_params(axis='x', labelsize=fs)
+    # #ax_sum.tick_params(axis='y', labelsize=fs)
+    # ax_sum.set_xscale('log')  # Logarithmische Skala für die Summen
+    # ax_sum.tick_params(axis='y', left=False, labelleft=False)
+    # ax_sum.set_xlim(1, 600)
+    
+    
+    AeqC = False
+    if AeqC:
+        a=2
+        ax_main.axhline(y=a, color='r', linestyle = '-', lw=1.5, label='A = '+str(a)+'C  |  C = '+str(a)+'A',)
+        ax_main.axhline(y=1/a, color='r', linestyle = '-', lw=1.5,)
+        a=2
+        ax_sum.axhline(y=a, color='r', linestyle = '-', lw=1.5,)
+        ax_sum.axhline(y=1/a, color='r', linestyle = '-', lw=1.5)
+        ax_main.legend(fontsize=fs, loc="upper left",ncol=2)
+    
+    # Achsenlimits setzen
+    limits = args.limits if args.limits else [x_edges[0], x_edges[-1], y_edges[0], y_edges[-1]]
+    ax_main.set_xlim(limits[0], limits[1])
+    ax_main.set_ylim(limits[2], limits[3])
+    #ax_sum.set_xlim(limits[0], limits[1])  # Gleiche x-Achse wie der Hauptplot
+    
+    # Layout optimieren
+    plt.tight_layout(rect=[0, 0, 0.95, 1])  # Platz für Colorbar rechts
+    save_or_show_plot([])
+
+
+###############################################################################    
+    
+def read_data(file, end):
+    return pd.read_csv(f"{file}.{end}", sep='\s+', skiprows=0)
+ 
+    
+###############################################################################
+
+def read_Itime(file, end):    
+    f = file + '.Itime'
+    Itime = None
+    while len(f) > 0:
+        if os.path.isfile(f+'.Itime'):
+            Itime = int(pd.read_csv(f+'.Itime', sep='\s+', skiprows=0, names = ['time'])['time'][0])
+            f=''
+        else:
+            f=f[:-1]
+    return Itime
+
+    
+###############################################################################
+    
+def save_or_show_plot(channels):
+    if args.export:
+        if args.export[0] != 'f':
+            name = args.export[0] 
+        else:
+            name = f'plots/{end}_{file.split("/")[-1]}'
+            if channels:
+                name += '_ch-' + ''.join(channels)
+        plt.tight_layout()
+        plt.savefig(f"{name}.{args.export[1] if len(args.export) > 1 else 'png'}")
+        print(f'saved as {name}')
+    plt.show()  
+   
+    
+################################
+def cutA1(A1):
+    arr=[]
+    for val in A1:
+        if val<8: arr.append(val+11)
+        else: arr.append((A1-8)/0.008)
+    return arr
+    
+def cutA2(A1):
+    arr=[]
+    for val in A1:
+        if val<8: arr.append((A1-7)/0.15)
+        else: arr.append(0.008*A1+11)
+    return arr
+###############################################################################
+
+def constant(x, c):
+    #return np.full_like(x,c)
+    return c
+
+###############################################################################
+
+if __name__ == '__main__':
+    main()

xplot.py

@@ -12,7 +12,7 @@ import os.path
 from matplotlib.colors import ListedColormap, Normalize
 from scipy.optimize import curve_fit
 from matplotlib.ticker import FuncFormatter
- 
+
 parser = argparse.ArgumentParser(description="This script plots histogramms of given data ending on '.2dhist'")
 parser.add_argument("-v", "--verbose", choices = ['0','1','2'], help = "when zero only emit serious error messages; when one also emit warnings, when two emit informational messages")
 parser.add_argument("filepath", type = str, help="path/file were data is stored in")
@@ -96,10 +96,7 @@ def do_1dhist():
     #limits = args.limits if args.limits else [-100, 10000, 0.1, 100000]
     limits = args.limits if args.limits else [0, 100, 10, 10000]
     
-    border = False
+    plt.axvline(x=2500, color='black', linestyle='--', linewidth=2, label='Gain switch (2500 $mV$)')
-    
-    if border:
-        plt.axvline(x=2500, color='black', linestyle='--', linewidth=2, label='Gain switching at 2500 $mV$')
         
     #plt.axvline(x=0.5, color='r', linestyle = '-', lw=1.5,)
     #plt.axvline(x=2, color='r', linestyle = '-', lw=1.5,)
@@ -250,7 +247,7 @@ def do_2dhist():
 
             #plt.plot([7.5, 10000, 10000, 7.5, 7.5], [7.5, 7.5, 10000, 10000, 7.5], color='r', linewidth=1, label = 'zoom')
             
-
+        plt.axhline(y=1, color='r', linestyle = '-', lw=1.5,)
             
         if ACcuts:
             # plt.axvline(x=7.5, color='r', linestyle = '-', lw=3, label = "thr(C1) = 7.5 mV")
@@ -305,10 +302,11 @@ def do_2dhist():
      
     else:
         #cb = plt.pcolormesh(x_edges, y_edges,c, cmap = cmap, shading='nearest', vmin=vmin, vmax = 2000,norm='log')  
-        cb = plt.pcolormesh(x_edges, y_edges,c, cmap = cmap, shading='auto', vmin=vmin, norm='log')
+        cb = plt.pcolormesh(x_edges, y_edges,c, cmap = cmap, shading='auto', vmin=vmin, vmax = 100, norm='log')
         #cb = plt.pcolormesh(np.multiply(x_edges,f), y_edges,c, cmap = cmap, shading='nearest', vmin=vmin, vmax=30, norm='log')    
         cbar = plt.colorbar(cb, orientation = 'vertical', shrink = 0.95, pad = 0.1)
-        
+      
+    
         if Lowcal:
             a = 40
             b = 130
@@ -412,7 +410,7 @@ def plot_layout(title, energy, limits, x, y, cbar):
         cbar.ax.tick_params(labelsize=fs)
     
     if plt.gca().get_legend_handles_labels()[0]:  # Gibt eine Liste von Künstlern (Handles) und Labels zurück
-        plt.legend(fontsize=fs-1, loc="upper right",ncol=1)
+        plt.legend(fontsize=fs-1, loc="upper right",ncol=4)
         #plt.legend(fontsize=fs+5, loc="upper right")