Update HV-ramp, add HV_hists

HV_hists.py

@@ -0,0 +1,123 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import glob
+import os
+import scienceplots 
+
+plt.style.use(['science'])
+fs = 26
+
+plt.rcParams["figure.figsize"] = (12,6)
+plt.rcParams['xtick.labelsize'] = fs-2
+plt.rcParams['ytick.labelsize'] = fs-2
+plt.rcParams['axes.labelsize'] = fs
+plt.rcParams['axes.titlesize'] = fs+2
+plt.rcParams['legend.fontsize'] = fs-6
+plt.rcParams['savefig.bbox'] = 'tight'
+
+# ------------------------------------------------------------
+# INSERT YOUR PRECOMPUTED HV VALUES HERE
+# ------------------------------------------------------------
+HV_labels_old = [1046, 1040, 1020, 998, 977, 956, 935, 913, 892, 870, 849]
+#HV_labels_old = [1046, 1040, 1020, 998, 977, 956, 935, 913, 892, 870]
+HV_labels = [970.2, 964.6, 945.6, 926.3, 906.6, 886.7, 866.9, 847.0, 827.1, 806.9, 786.7]
+HV_labels = [970.2, 964.6, 945.6, 926.3, 906.6, 886.7, 866.9, 847.0, 827.1, 806.9, 786.7]
+
+# ------------------------------------------------------------
+# Histogram folder
+# ------------------------------------------------------------
+hist_path = "ramp_hists/"
+
+# ------------------------------------------------------------
+# Helper to read histogram file
+# ------------------------------------------------------------
+def read_histogram(file):
+    data = np.loadtxt(file, comments='#', skiprows=1)
+    x = data[:, 0]
+    y = data[:, 1]
+    return x, y
+
+def read_itime(file):
+    """Lese Messzeit aus .Itime Datei (in s)"""
+    with open(file) as f:
+        t_total = float(f.read().strip())
+    return t_total
+
+# ------------------------------------------------------------
+# Collect files
+# ------------------------------------------------------------
+# files_1d  = sorted(glob.glob(hist_path + "*_mV_ch-B.1dhist"))
+# files_old = sorted(glob.glob(hist_path + "*_Bhist.hist"))
+files_1d = []
+files_old = []
+files_itime = []
+
+for hv in HV_labels_old:
+    # Suche die passende 1dhist-Datei für diese HV-Stufe
+    match_1d = glob.glob(os.path.join(hist_path, f"*{hv}V*_mV_ch-B.1dhist"))
+    if match_1d:
+        files_1d.append(match_1d[0])
+    else:
+        print(f"Warning: No 1dhist file found for {hv} V")
+
+    # Suche die passende Bhist-Datei
+    match_old = glob.glob(os.path.join(hist_path, f"*{hv}V*_Bhist.hist"))
+    if match_old:
+        files_old.append(match_old[0])
+    else:
+        print(f"Warning: No Bhist file found for {hv} V")
+
+    # Itime
+    match_time = glob.glob(os.path.join(hist_path, f"*{hv}V*.Itime"))
+    if match_time:
+        files_itime.append(match_time[0])
+    else:
+        print(f"Warning: No Itime file found for {hv} V")
+        files_itime.append(None)  # Platzhalter
+
+# ------------------------------------------------------------
+# Plot 1 — mV_ch-B.1dhist
+# ------------------------------------------------------------
+plt.figure()
+
+for file, hv ,itime_file in zip(files_1d, HV_labels,files_itime):
+    x, y = read_histogram(file)
+    bin_width = x[1] - x[0]
+
+    if itime_file is not None:
+        t_total_s = read_itime(itime_file)
+        t_total_h = t_total_s / 3600
+        y = y / (t_total_h * bin_width)  # counts / h / mV
+    else:
+        print(f"Warning: Using raw counts for {hv} V (no Itime)")
+
+    plt.step(x, y, where='mid', label=f"{hv:.0f} V")
+
+plt.yscale("log")
+plt.xlabel("Signal [mV]")
+plt.xlim(-10,1000)
+plt.ylabel("Counts [1/h/mV]")
+plt.title("Ramp Histograms (mV channel B)")
+plt.grid(True)
+plt.legend(frameon=True, framealpha=0.85)
+plt.tight_layout()
+
+# ------------------------------------------------------------
+# Plot 2 — Bhist.hist
+# ------------------------------------------------------------
+plt.figure()
+
+for file, hv in zip(files_old, HV_labels):
+    x, y = read_histogram(file)
+    plt.step(x, y, where='mid', label=f"{hv:.0f} V")
+
+plt.yscale("log")
+plt.xlabel("Signal [mV]")
+plt.xlim(-10,1000)
+plt.ylabel("Counts [1/h/mV]")
+plt.title("Ramp Histograms (Bhist)")
+plt.grid(True)
+plt.legend(frameon=True, framealpha=0.85)
+plt.tight_layout()
+
+plt.show()

HV_photoelectrons.py

@@ -0,0 +1,97 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.optimize import curve_fit
+from matplotlib.ticker import ScalarFormatter
+from scipy.constants import physical_constants
+
+# Elementarladung
+e = physical_constants['elementary charge'][0]
+
+# Plot-Stil
+plt.rcParams["figure.figsize"] = (15,8)
+plt.rcParams['xtick.labelsize'] = 15
+plt.rcParams['ytick.labelsize'] = 15
+plt.rcParams['axes.labelsize'] = 20
+plt.rcParams['axes.titlesize']=25
+plt.rcParams['legend.fontsize']=20
+plt.rcParams['savefig.bbox']='tight'
+
+# --- Datenklasse (wie im ersten Skript) ---
+class data():
+    def __init__(self, path='test.EI'):
+        self.path = path[:-3]
+        self.mV = 14000
+        self.Itime = 0
+        self.ei_lines = 0
+        self.list = [[] for _ in range(18)]
+        self.keylist = ['A1L','A1H','E1','A2L','A2H','E2','C1L','C1H','E3','C2L','C2H','D1L','EL','EH','D1H','D2L','D2H','B']
+        self.data = {key: None for key in self.keylist}
+
+        with open('data/EIs/' + path, 'r') as f:
+            writing = False
+            for line in f:
+                parts = line.split()
+                if len(parts) == 0:
+                    continue
+                if parts[0] == 'HC':
+                    writing = not writing
+                    if writing == False:
+                        self.Itime += 10
+                        self.ei_lines = 0
+                if parts[0] == 'EI' and parts[-1] != 'EOF' and len(parts) == 59 and writing and int(parts[4].split(':')[0]) > 10000:
+                    self.ei_lines += 1
+                    for i in range(18):
+                        self.list[i].append(int(parts[5 + 3*i])/self.mV)
+        
+        for i in range(18):
+            if self.ei_lines != 0:
+                self.list[i] = np.array(self.list[i][:-self.ei_lines])
+            else:
+                self.list[i] = np.array(self.list[i])
+
+    # Histogramm
+    def hist(self, resV=20):
+        minV, maxV = -10, 3600
+        bins = int((maxV - minV)/resV)
+        for i, key in enumerate(self.data):
+            hist, X = np.histogram(self.list[i], bins, (minV, maxV))
+            hist = hist / (resV * (self.Itime/3600))
+            self.data[key] = hist, X[:-1]
+        # noM
+        mask = (self.list[14] < 50) & (self.list[16] < 50)
+        hist, X = np.histogram(self.list[17][mask], bins, (minV, maxV))
+        hist = hist / (resV * (self.Itime/3600))
+        self.data['noM'] = hist, X[:-1]
+
+    # Gauss-Funktion
+    def gaus_function(self, x, a, x0, sigma):
+        return a*(1/(np.sqrt(2*np.pi)*sigma))*np.exp(-(x-x0)**2/(2*sigma**2))
+
+    # Fit auf Gauss
+    def fit_gauss(self, key='noM', range_mask=None, startpoints=[500,80,20]):
+        x = self.data[key][1]
+        y = self.data[key][0]
+        if range_mask is not None:
+            x = x[range_mask]
+            y = y[range_mask]
+        self.data['popt_gauss'], self.data['pcov_gauss'] = curve_fit(self.gaus_function, x, y, p0=startpoints)
+        self.data['fitX'] = np.arange(min(x), max(x), 1)
+        self.data['fitY'] = self.gaus_function(self.data['fitX'], *self.data['popt_gauss'])
+        print(f"Gauss Fit: mu={self.data['popt_gauss'][1]:.2f}, sigma={self.data['popt_gauss'][2]:.2f}, sigma/mu={self.data['popt_gauss'][2]/self.data['popt_gauss'][1]:.3f}")
+
+# --- Chaos Flight einlesen ---
+chaos = data('chaos_flight.EI')
+chaos.hist()
+chaos.fit_gauss()
+
+# --- Plot: Counts/h/mV vs Photoelectrons ---
+plt.figure()
+plt.step(chaos.data['noM'][1], chaos.data['noM'][0], label='Data')
+plt.plot(chaos.data['fitX'], chaos.data['fitY'], 'r--', label='Gauss Fit')
+plt.xlabel('Photoelectrons (PE)')
+plt.ylabel('Counts [1/mV/h]')
+plt.yscale('log')
+plt.grid()
+plt.legend()
+plt.title('CHAOS Flight Histogram')
+plt.show()

HV_ramp.py

@@ -0,0 +1,177 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import scienceplots 
+
+plt.style.use(['science'])
+fs = 26
+
+# ------------------------------------------------------------
+# Plot style
+# ------------------------------------------------------------
+plt.rcParams["figure.figsize"] = (12,6)
+plt.rcParams['xtick.labelsize'] = fs-2
+plt.rcParams['ytick.labelsize'] = fs-2
+plt.rcParams['axes.labelsize'] = fs
+plt.rcParams['axes.titlesize'] = fs+2
+plt.rcParams['legend.fontsize'] = fs-4
+plt.rcParams['savefig.bbox'] = 'tight'
+
+# ------------------------------------------------------------
+# Files
+# ------------------------------------------------------------
+file_ramp   = 'data_justus/2024-12-28-chaos-HVramp-chk-25mV-3_25h.EI'
+file_flight = 'chaos_sd_float_new.EI'
+
+# Split ramp files (already time-cut where needed)
+ramp_files = [
+    'CHAOS_ramp_1046V.EI',
+    'CHAOS_ramp_1040V.EI',
+    'CHAOS_ramp_1020V.EI',
+    'CHAOS_ramp_998V.EI',
+    'CHAOS_ramp_977V.EI',
+    'CHAOS_ramp_956V.EI',
+    'CHAOS_ramp_935V.EI',
+    'CHAOS_ramp_913V.EI',
+    'CHAOS_ramp_892V.EI',
+    'CHAOS_ramp_870V.EI',
+    'CHAOS_ramp_849V_5h.EI'
+]
+
+ramp_path = "data_justus/"
+
+# ------------------------------------------------------------
+# HV conversion constants
+# ------------------------------------------------------------
+ADC_MID   = 0x8000
+ADC_RANGE = 5.12
+HV_GAIN   = 281
+
+# ------------------------------------------------------------
+# Read HC data
+# ------------------------------------------------------------
+def read_HC(file):
+    t, hv_adc = [], []
+    with open(file) as f:
+        for line in f:
+            p = line.split()
+            if p and p[0] == 'HC':
+                t.append(float(p[1]))
+                hv_adc.append(float(p[11]))
+
+    t = np.array(t)
+    hv_adc = np.array(hv_adc)
+
+    hv = (hv_adc / ADC_MID - 1) * ADC_RANGE * HV_GAIN
+    time_h = (t - t[0]) / 3600
+
+    return time_h, hv
+
+# ------------------------------------------------------------
+# Find HV steps (common function for both plots)
+# ------------------------------------------------------------
+def find_HV_steps(time, HV, threshold=5):
+    steps = []
+    current = [0]
+
+    for i in range(1, len(HV)):
+        if abs(HV[i] - HV[i-1]) > threshold:
+            steps.append(current)
+            current = [i]
+        else:
+            current.append(i)
+
+    steps.append(current)
+    HV_means = [np.mean(HV[s]) for s in steps]
+
+    return steps, HV_means
+
+# ------------------------------------------------------------
+# Plot 1: Original HV Ramp + Flight
+# ------------------------------------------------------------
+time, HV = read_HC(file_ramp)
+time_f, HV_f = read_HC(file_flight)
+HV_f_mean = np.mean(HV_f)
+
+plt.figure()
+
+plt.plot(time, HV, label='HV Ramp', color='blue')
+plt.plot(time_f, HV_f, color='tab:orange', label=f'Float Data; HV mean = {HV_f_mean:.0f} V')
+
+steps, HV_steps = find_HV_steps(time, HV)
+# print("HV_labels = [")
+# for m in HV_steps:
+#     print(f"    {m:.1f},")
+# print("]")
+
+for s, m in zip(steps, HV_steps):
+    t0, t1 = time[s[0]], time[s[-1]]
+    plt.hlines(m, t0, t1, colors='k', linestyles='--')
+    plt.text((t0+t1)/2, m+5, f"{m:.0f} V", ha='center', fontsize=fs-8)
+
+plt.hlines(HV_f_mean,
+           time[0],
+           time[-1],
+           colors='tab:orange',
+           linestyles='--')
+
+plt.xlabel('Time since start [h]')
+plt.ylabel('High Voltage [V]')
+plt.title('HV Ramp + Float')
+plt.xlim(0, 25)
+plt.ylim(None, 1000)
+plt.margins(x=0)
+plt.grid(True)
+plt.legend(frameon=True, framealpha=0.85)
+plt.savefig("Bplots/HV_ramp.pdf")
+plt.savefig("Bplots/HV_ramp.png")
+plt.tight_layout()
+
+
+# # ------------------------------------------------------------
+# # Plot 2: Reconstructed Ramp from Split Files
+# # ------------------------------------------------------------
+# plt.figure()
+
+# time_offset = 0
+# all_times = []
+# all_HV = []
+
+# for f in ramp_files:
+#     t, hv = read_HC(ramp_path + f)
+
+#     t_shifted = t + time_offset
+#     all_times.append(t_shifted)
+#     all_HV.append(hv)
+
+#     time_offset = t_shifted[-1]
+
+# time_concat = np.concatenate(all_times)
+# HV_concat   = np.concatenate(all_HV)
+
+# plt.plot(time_concat, HV_concat, label='HV Ramp', color='blue')
+# plt.plot(time_f, HV_f, color='tab:orange', label='Flight Data')
+
+# steps_rec, HV_steps_rec = find_HV_steps(time_concat, HV_concat)
+
+# for s, m in zip(steps_rec, HV_steps_rec):
+#     t0, t1 = time_concat[s[0]], time_concat[s[-1]]
+#     plt.hlines(m, t0, t1, colors='k', linestyles='--')
+#     plt.text((t0+t1)/2, m+5, f"{m:.0f} V", ha='center')
+
+# plt.hlines(HV_f_mean,
+#            time_concat[0],
+#            time_concat[-1],
+#            colors='tab:orange',
+#            linestyles='--',
+#            label=f'Flight HV mean = {HV_f_mean:.0f} V')
+
+# plt.xlabel('Time since start [h]')
+# plt.ylabel('High Voltage [V]')
+# plt.title('Reconstructed HV Ramp from Split Files')
+# plt.xlim(left=0)
+# plt.margins(x=0)
+# plt.grid(True)
+# plt.legend()
+# plt.tight_layout()
+
+plt.show()