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,7 +102,6 @@ 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
@@ -119,13 +118,13 @@ def get_values(l,T):
 
 deltacut = 0
 
-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','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=[]
@@ -145,8 +144,8 @@ folder = "Bhists"
 
 minX = 2
 maxX = 10000
-minY = 0.1
-maxY = 2000
+minY = 0.02
+maxY = 20
 resV = 50*80*resV
 resX = resV/0.00362/polynom(15,*pol)
 
@@ -161,7 +160,6 @@ 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

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

@@ -1,41 +0,0 @@
-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,7 +191,6 @@ def do_2dhist():
     Ecuts = False
     Lowcal = False
     mask_He = False
-    Bcuts = True 
     
     if end == '2dhist':
         vmin = 1
@@ -290,10 +289,7 @@ def do_2dhist():
             plt.fill_between(xx[a:], 1/0.01*xx[a:],100000, color='b', alpha=alpha)
             
             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
@@ -306,7 +302,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 = 1000,norm='log')  
+        #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, 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')