add two skripts from justus

CHAOSa.py

@@ -603,15 +603,15 @@ def generate_2d_hist_logxy():
     #maxY=200000
     
     #A1 vs A2
-    maxX = 100000
-    maxY = 100000
+    #maxX = 100000
+   # maxY = 100000
     #maxY = 10
     
     #maxX=3500
     #maxY=3500
     
-    minY = 0.01
-    maxY = 100
+    #minY = 0.01
+    #maxY = 100
 
     bin_x = int((maxX-minX)/resX)
     bin_y = bin_x+1

CHAOScustom.py

@@ -102,6 +102,7 @@ def get_values_fish(l,T):
 def get_values_other(l,T):
     y = min(value_HL(l,"A",u),value_HL(l,"C",u))/1000
     #y = value_HL(l,"A",u)/1000
+    y = value_HL(l,"B",u)
     x = value_HL(l,"D1",u)+value_HL(l,"D2",u)
     x2 = get_MeV(x,T)
     return x2,y
@@ -118,13 +119,13 @@ def get_values(l,T):
 
 deltacut = 0
 
-allcuts = [['A','C','AeqC','D'],['A','B','C','AeqC','D'],['A','nB','C','AeqC','D'],['A','rB','C','AeqC','D']]
+allcuts = [[],['D'],['A','C','D'],['A','C','AeqC','D']]#,['A','B','C','AeqC','D'],['A','nB','C','AeqC','D'],['A','rB','C','AeqC','D']]
 #allcuts = [['A','Bsim','C','AeqC','D'],['A','nBsim','C','AeqC','D']]
 #allcuts = [['A','D'],['A','B','D'],['A','nB','D']]
-allcuts = [['A','B','C','AeqC','D'],['A','nB','C','AeqC','D'],['A','C','AeqC','D']]
+#allcuts = [['A','B','C','AeqC','D'],['A','nB','C','AeqC','D'],['A','C','AeqC','D']]
 #allcuts = [['A','C','AeqC','D','E123'],['A','C','AeqC','D','E']]
 #allcuts = [[],['A','C'],['D'],['A','C','D'],['A','C','AeqC','D']]
-allcuts=[['A','B','C','AeqC','D'],['A','nB','C','AeqC','D'],['A','rB','C','AeqC','D'],['A','B1','C','AeqC','D'],['A','B2','C','AeqC','D'],['A','B3','C','AeqC','D']]
+#allcuts=[['A','B','C','AeqC','D'],['A','nB','C','AeqC','D'],['A','rB','C','AeqC','D'],['A','B1','C','AeqC','D'],['A','B2','C','AeqC','D'],['A','B3','C','AeqC','D']]
 
 #trig = ['D1H','D2H']
 trig=[]
@@ -144,8 +145,8 @@ folder = "Bhists"
 
 minX = 2
 maxX = 10000
-minY = 0.02
-maxY = 20
+minY = 0.1
+maxY = 2000
 resV = 50*80*resV
 resX = resV/0.00362/polynom(15,*pol)
 
@@ -160,6 +161,7 @@ if allother:
         #filename = "chaos_Cut"+ "".join(cut)+"_trig-"+ "".join(trig)+"~min(A,C)-D.2dhist2log"
         #filename = "fullsimProton_Cut"+ "".join(cut)+"~min(A,C)-D.2dhist2log"
         filename = "chaosFloat_Cut"+ "".join(cut)+"~min(A,C)-D.2dhist2log"
+        filename = "chaosFloat_Cut"+ "".join(cut)+"~B-D.2dhist2log"
         do_hist(cut,trigchans,'other') 
         
 if allfish:

justus_data_eval_CHAOS.py

@@ -0,0 +1,254 @@
+import numpy as np
+import pickle
+import matplotlib.pyplot as plt
+from scipy.optimize import curve_fit
+from matplotlib.ticker import ScalarFormatter
+from scipy.constants import c, physical_constants
+e = physical_constants['elementary charge'][0]
+
+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'
+
+
+class data():
+
+    def __init__(self, path='test.EI'):
+
+        self.path = path[:-3]
+        self.mV = 14000
+        self.ei_lines = 0
+        self.Itime = 0
+        self.list = [[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]]
+        #Die Belegungen der Kanälen sind in Pierres Dokumentation gut dargestellt.
+        self.keylist = ['A1L','A1H','E1','A2L','A2H','E2','C1L','C1H','E3','C2L','C2H','D1L','EL','EH','D1H','D2L','D2H','B']
+        self.data ={}
+
+        #Erzeugt Dictonary mit den ausgelesenen Kanälen:
+        for key in self.keylist:
+            self.data[key] = None
+
+        if self.list[0] == []:
+            with open('data/EIs/'+path,'r') as f:
+                writing = False
+
+                for line in f:
+                    if line.split()[0]=='HC' and writing == False:
+                        writing = True
+                    elif line.split()[0]=='HC':
+                        self.Itime += 10
+                        self.ei_lines = 0
+                    if line.split()[0]=='EI'and line.split()[-1]!='EOF' and len(line.split())==59 and writing == True  and int((line.split()[4]).split(':')[0])>10000:
+                        self.ei_lines += 1
+                        self.list[0].append(int(line.split()[5])/self.mV)
+                        self.list[1].append(int(line.split()[8])/self.mV)
+                        self.list[2].append(int(line.split()[11])/self.mV)
+                        self.list[3].append(int(line.split()[14])/self.mV)
+                        self.list[4].append(int(line.split()[17])/self.mV)
+                        self.list[5].append(int(line.split()[20])/self.mV)
+                        self.list[6].append(int(line.split()[23])/self.mV)
+                        self.list[7].append(int(line.split()[26])/self.mV)
+                        self.list[8].append(int(line.split()[29])/self.mV)
+                        self.list[9].append(int(line.split()[32])/self.mV)
+                        self.list[10].append(int(line.split()[35])/self.mV)
+                        self.list[11].append(int(line.split()[38])/self.mV)
+                        self.list[12].append(int(line.split()[41])/self.mV)
+                        self.list[13].append(int(line.split()[44])/self.mV)
+                        self.list[14].append(int(line.split()[47])/self.mV)
+                        self.list[15].append(int(line.split()[50])/self.mV)
+                        self.list[16].append(int(line.split()[53])/self.mV)
+                        self.list[17].append(int(line.split()[56])/self.mV)
+                for i, list  in enumerate(self.list):
+                    if self.ei_lines != 0:
+                        self.list[i] = list[:-self.ei_lines]
+                    self.list[i] = np.array(self.list[i])
+                f.close()
+
+        
+    def hist(self, resV=20):
+        resV = resV
+        minV = -10
+        maxV = 3600
+        bins=int((maxV-minV)/resV)
+        for i, key in enumerate(self.data):
+            if i < len(self.list):
+                hist, X = np.histogram(self.list[i], bins, (minV, maxV))
+                hist = hist/(resV*(self.Itime/3600))
+                self.data[key] = hist, X[:-1]
+
+        # Maske für wenn BGO ein Myon gesehen hat. Dafür sollen entweder D1 und D2 High gain beide 65mV gesehen haben. die low gain kanäle sind zu dicht am rauschen.
+        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]
+
+        #Cut für Myon in BGO und in Detektor Stufe A und C gesehen:
+        hist, X = np.histogram(self.list[17][self.cut()], bins, (minV, maxV))
+        hist = hist/(resV*(self.Itime/3600))
+        self.data['muon'] = hist, X[:-1]
+
+
+
+    def cut(self, muon = True):
+        #BGO:
+        CUT_BGO = (self.list[14]>50) & (self.list[16]>50) & (self.list[16]>self.list[14]*0.8-40) & (self.list[16]<self.list[14]*1.2+40)
+        #SSDA:
+        CUT_SSDA = (self.list[1]>13) | (self.list[4]>17)
+        #SSDB
+        CUT_SSDC = (self.list[7]>13) | (self.list[10]>17)
+        return (CUT_BGO & CUT_SSDA & CUT_SSDC)
+
+    def plot(self, plots=['B']):
+        for key in plots:
+            plt.step(self.data[key][1], self.data[key][0], label=key)
+        try:
+            plt.plot(self.data['fitX'], self.data['fitY'], label='fit', linestyle='--',color='black')
+            plt.plot(np.arange(-10,3500,1), self.data['fit_gauss'], label='fit_gauss', color='red')
+        except:
+            pass
+        plt.yscale('log')
+        plt.ylim(1,2*10**4)
+        plt.xlim(-10,3500)
+        plt.grid()
+        plt.xlabel('Signal [mV]')
+        plt.ylabel('Zählrate [1/mV/h]')
+        plt.legend(loc='upper right')
+        plt.title(self.path)
+        plt.show()
+
+    def fit_function(self,x, a, x0, sigma ,b,c,b1,c1):
+        #guss_pedastal= a1*np.exp(-(x-x1)**2/(2*sigma1**2))# a1, x1, sigma1
+        return  (a*(1/(np.sqrt(2*np.pi)*sigma))*np.exp(-(x-x0)**2/(2*sigma**2))+ b*np.exp(-c*x)+b1*np.exp(-c1*x))
+    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)))
+    
+    def fit(self, range=None, startpoints=[1*10**4,80,200,500,0.02,5,0.0005]):
+        if range == None:
+            range=(self.data['noM'][1]>=25)&(self.data['noM'][1]<700)        
+        self.data['popt'], self.data['pcov'] = curve_fit(self.fit_function,self.data['noM'][1][range], self.data['noM'][0][range], p0=startpoints)
+        self.data['fitX'] = np.arange(24,3100,1)
+        self.data['fitY'] = self.fit_function(self.data['fitX'], *self.data['popt'])
+        self.data['fit_gauss'] = self.gaus_function(np.arange(-10,3500,1), *self.data['popt'][:3])
+        print('relative Breite: ', self.data['popt'][2]/self.data['popt'][1])
+
+
+try:
+    with open('processed_data_HV_RAMP','rb') as file:
+        loaded_data = pickle.load(file)
+        print('Vorhandene Daten wurden eingelesen')
+except:
+    print('Noch keine prozessierten Daten vorhanden')
+    
+
+'''
+In plot2 werden Plots verschiedener Objekte (also verschiedener Datensätze/Messungen) erzeugt.
+Man übergibt mit plots eine Liste der gewünschten Plots. Pro Plot steht in der Liste ein Tupel.
+Im ersten Eintrag des Tupels steht das Histogramm, im zweiten das label. Fügt man dem Plot einen Tite hinzu, so wird er in Plots/ abgeschpeichert.
+'''
+
+def plot_leer(plots=[], labels=['PMT', 'PMT', 'PMT'], title=None, colors=['tab:blue', 'tab:orange', 'tab:red']):
+        for index, key in enumerate(plots):
+            plt.step(key.data['noM'][1], key.data['noM'][0], label=labels[index])
+            plt.plot(key.data['fitX'], key.data['fitY'], linestyle='--',color='black')
+            #plt.plot(np.arange(-10,3500,1), key.data['fit_gauss'], label= r'Fit Gauss: $\mu = $'+str(int(key.data['popt'][1]))+r'; $\sigma = $'+str(int(key.data['popt'][2]))+r'; $\frac{\sigma}{\mu} = $'+ str(np.round(key.data['popt'][2]/key.data['popt'][1],3)), color=colors[index])
+        plt.yscale('log')
+        plt.ylim(55,5*10**3)
+        plt.xlim(10,150)
+        plt.grid()
+        plt.xlabel('Signal [mV]')
+        plt.ylabel('Counts [1/mV/h]')
+        plt.legend(loc='upper right',ncol=2)
+        if title != None:
+            #plt.title(title)
+            plt.savefig('Plots/CHAOS/'+title,dpi=1000)
+        plt.show()
+
+#Umrechnungsfaktor: 5.12/0x8000*303
+#Cut am BGO: summe >50 jeder >30
+
+'''HV_1046 = data('CHAOS_ramp_1046V.EI')
+HV_1040 = data('CHAOS_ramp_1040V.EI')
+HV_1020 = data('CHAOS_ramp_1020V.EI')
+HV_998 = data('CHAOS_ramp_998V.EI')
+HV_977 = data('CHAOS_ramp_977V.EI')
+HV_956 = data('CHAOS_ramp_956V.EI')
+HV_935 = data('CHAOS_ramp_935V.EI')
+HV_913 = data('CHAOS_ramp_913V.EI')
+HV_892 = data('CHAOS_ramp_892V.EI')
+HV_872 = data('CHAOS_ramp_870V.EI')
+HV_849 = data('CHAOS_ramp_849V.EI')'''
+
+#plot_leer(plots=[HV_1046, HV_1040, HV_1020, HV_998, HV_977, HV_956, HV_935, HV_913, HV_892, HV_872, HV_849], labels=['1046 V', '1040 V', '1020 V', '998 V', '997 V', '956 V', '935 V', '913 V', '892 V', '872 V', '849 V'])
+'''
+
+
+#Verstärkungplot:
+def plot_gain():
+    X = []
+    Y = []
+    for key in loaded_data:
+        X.append(key)
+        Y.append(loaded_data[key].data['popt'][1])
+    Y = np.array(Y)
+    print(Y)
+    gain = Y/(1000*15*10**10*e)
+    plt.errorbar(X,Y,yerr=0.1*Y, capsize=3, fmt='+')    
+    plt.xlabel('PMT Voltage [V]')
+    plt.ylabel('Signal SE [mV]')
+    plt.text(844,32,'CHAOS \n Flight', fontsize = 15, bbox=dict(boxstyle="round", edgecolor="black", facecolor="lightblue"))
+    plt.yscale('log')
+    X[1]=1035
+    plt.xticks(X)
+    plt.yticks([28.0,30, 40, 50, 60])
+    plt.gca().yaxis.set_major_formatter(ScalarFormatter())  # Nutzt Dezimalformat
+    plt.ticklabel_format(style='plain', axis='y')  # Unterdrückt wissenschaftliches Format
+    plt.grid()
+    plt.savefig('Plots/CHAOS/'+'gain_in_MV',dpi=1000)
+    plt.show()'''
+# Verstärkungplot:
+def plot_gain():
+    X = []
+    Y = []
+    for key in loaded_data:
+        X.append(key)
+        Y.append(loaded_data[key].data['popt'][1])
+    Y = np.array(Y)
+    print(Y)
+    gain = Y / (1000 * 15 * 10**10 * e)
+    v_datasheet = np.array([850,950,1020])
+    gain_datasheet = np.array([6*10**5,9*10**5,2.2*10**6])*(1000 * 15 * 10**10 * e)
+    
+    fig, ax1 = plt.subplots()
+
+    # Linke Achse (Signal SE [mV])
+    ax1.errorbar(X, Y, yerr=0.1 * Y, capsize=3, fmt='+')
+    ax1.errorbar(v_datasheet,gain_datasheet, capsize=3, fmt='+')
+    ax1.set_xlabel('PMT Voltage [V]')
+    ax1.set_ylabel('Signal SE [mV]')
+    ax1.set_yscale('log')
+    ax1.set_xticks(X)
+    ax1.set_yticks([28.0, 30, 40, 50, 60])  # Beispielwerte, anpassen, falls nötig
+    ax1.yaxis.set_major_formatter(ScalarFormatter())
+    ax1.ticklabel_format(style='plain', axis='y')
+    ax1.grid()
+
+    # Text hinzufügen
+    ax1.text(844, 32, 'CHAOS \n Flight', fontsize=15, bbox=dict(boxstyle="round", edgecolor="black", facecolor="lightblue"))
+
+    # Rechte Achse (Gain)
+    ax2 = ax1.twinx()
+    ax2.set_ylabel('Gain')
+    ax2.set_yscale('log')
+    ax2.set_ylim(ax1.get_ylim())  # Gleiche Skalierung wie die linke Achse, damit Werte korrespondieren
+    gain_ticks = ax1.get_yticks()  # Tick-Werte der linken Achse holen
+    ax2.set_yticks(gain_ticks)
+    ax2.set_yticklabels([f"{tick / (1000 * 15 * 10**10 * e):.2e}" for tick in gain_ticks])  # Werte für Gain umrechnen und anzeigen
+
+    # Speichern und Plot anzeigen
+    plt.savefig('Plots/CHAOS/'+'gain_in_MV', dpi=1000)
+    plt.show()
+

justus_housekeepingCHAOS.py

@@ -0,0 +1,41 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+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'
+
+HV = np.load('HV.npy')
+time = np.load('time.npy')
+
+'''Zum erzeugen der Arrays wurde folgender Shell Kommand ausgeführt:
+HV = !grep HC data/EIs/2024-12-28-chaos-HVramp-chk-25mV-3.EI | awk '{print $12}'
+time = !grep  HC data/EIs/2024-12-28-chaos-HVramp-chk-25mV-3.EI | awk '{print $2}'
+'''
+#Erzeugt aus Strings floats:
+HV = HV.astype(float)
+time = time.astype(float)
+
+#Rechnet ADC Werte in HV um:
+HV = (HV-float(0x8000))*(5.12/float(0x8000))*303
+#Macht aus falscher UNIX Zeit, Stunden seit Messbeginn:
+time = (time-time[0])/3600
+
+#Schmeißt alle fehlerhaften werte raus und alle die länger als 3 Stunden mit lettem Spannungswert gemessen haben:
+mask = (abs(HV) > 700) & (abs(HV) < 1400) & (time <=23)
+time = time[mask]
+HV = HV[mask]
+
+#Plotten:
+plt.plot(time,HV)
+plt.xlabel('Zeit seit Messbeginn [h]')
+plt.ylabel('Hochspannung [V]')
+plt.grid()
+plt.savefig('Plots/CHAOS/HV-ramp',dpi=600)
+plt.show()
+
+

xplot.py

@@ -191,6 +191,7 @@ def do_2dhist():
     Ecuts = False
     Lowcal = False
     mask_He = False
+    Bcuts = True 
     
     if end == '2dhist':
         vmin = 1
@@ -290,6 +291,9 @@ def do_2dhist():
             
             plt.fill_between(xx[a:], 1/0.01*xx[a:],np.maximum(0.01 * xx[a:], 12), color='orange', alpha=alpha, label = 'dual hits')
 
+        if Bcuts:
+               plt.axhline(y=17, color='b', linestyle = '-', lw=1.5) 
+               plt.axhline(y=43, color='b', linestyle = '-', lw=1.5)
             
         if mask_He:
             x1,y1= 64,0.6
@@ -302,7 +306,7 @@ def do_2dhist():
         cbar = plt.colorbar(contour, orientation = 'vertical', shrink = 0.95, pad = 0.1)
      
     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='nearest', vmin=vmin, vmax = 1000,norm='log')  
         cb = plt.pcolormesh(x_edges, y_edges,c, cmap = cmap, shading='auto', vmin=vmin, vmax = 100, norm='log')
         print(vmin)
         #cb = plt.pcolormesh(np.multiply(x_edges,f), y_edges,c, cmap = cmap, shading='nearest', vmin=vmin, vmax=30, norm='log')