Add code

calculate_edge_from_e.py

@@ -0,0 +1,27 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+from scipy.constants import m_e,c, elementary_charge
+
+
+class calculate_edge:
+    def __init__(self) -> None:
+        data = pd.read_csv('fitting_values_all.csv')
+        # remove column u_edge
+        data = data.drop(columns=['u_edge'])
+        self.energy_first_compton_edge = 569698
+        e = data['e'].to_numpy()
+        u_edge = self.calculate_calibration_first_edge(e)
+        print(np.round(u_edge,5))
+        # add to dataframe between the columns e and s 
+        data.insert(4, 'u_edge', np.round(u_edge,5))
+        # save now the dataframe
+        data.to_csv('fitting_values_all.csv', index=False)
+
+
+
+    def calculate_calibration_first_edge(self,e):
+        energie = self.energy_first_compton_edge* elementary_charge * (1 - 1 / (1 + 2 * self.energy_first_compton_edge * elementary_charge / (m_e * c**2)))
+        u = e/energie*1e3* elementary_charge
+        return u
+calculate_edge()

check_random_hits.py

@@ -0,0 +1,233 @@
+import numpy as np
+import pandas as pd
+import os
+import matplotlib.pyplot as plt
+import tqdm
+plt.rcParams.update({'font.size': 18})
+
+'''
+This script outputs the muon 1d histograms for all different combinations for my thesis setup.
+
+
+Kanalzuordnung:
+        TO1 = 0; thr[TO1] = 12; ch[0] = TO1; name[0]="TO1"
+        TU1 = 1; thr[TU1] = 12; ch[1] = TU1; name[1]="TU1"
+        TU2 = 2; thr[TU2] = 12; ch[2] = TU2; name[2]="TU2"
+        TO2 = 3; thr[TO2] = 12; ch[3] = TO2; name[3]="TO2"
+        TU3 = 4; thr[TU3] = 12; ch[4] = TU3; name[4]="TU3"
+        TO3 = 5; thr[TO3] = 12; ch[5] = TO3; name[5]="TO3"
+        TO4 = 6; thr[TO4] = 12; ch[6] = TO4; name[6]="TO4"
+        TU4 = 7; thr[TU4] = 12; ch[7] = TU4; name[7]="TU4"
+        TU5 = 8; thr[TU5] = 12; ch[8] = TU5; name[8]="TU5"
+        TO5 = 9; thr[TO5] = 12; ch[9] = TO5; name[9]="TO5"
+        TU7 = 10; thr[TU7] = 12; ch[10] = TU7; name[10] = "TU7"
+        TO6 = 11; thr[11] = 12; ch[11] = TO6; name[11] = "TO6"
+        TO7 = 12; thr[12] = 12; ch[12] = TO7; name[12] = "TO7"
+        TU6 = 13; thr[13] = 12; ch[13] = 13; name[13] = "TU6"
+        TU8 = 14; thr[14] = 12; ch[14] = 14; name[14] = "TU8"
+        TO8 = 15; thr[15] =12; ch[15] = 15; name[15] = "TO8"
+        B2 = 16; thr[16] =12; ch[16] = 16; name[16] = "B2"
+        B1 = 17; thr[17] = 12; ch[17] = 17; name[17] = "B1" 
+
+neue Kanalzuordnung:
+        oben:
+        TO1 = 3; thr[TO1] = 12; ch[0] = TO1; name[0]="TO1" 
+        TO2 = 5; thr[TO2] = 12; ch[1] = TO2; name[1]="TO2"
+        TO3 = 6; thr[TO3] = 12; ch[2] = TO3; name[2]="TO3"
+        TO4 = 9; thr[TO4] = 12; ch[3] = TO4; name[3]="TO4"
+        TO5 = 11; thr[TO5] = 12; ch[4] = TO5; name[4]="TO5"
+        TO6 = 15; thr[15] = 12; ch[5] = TO6; name[5] = "TO6"
+        TO7 = 0; thr[0] = 12; ch[6] = TO7; name[6] = "TO7"
+        TO8 = 12; thr[12] =12; ch[7] = 12; name[7] = "TO8"
+        unten:
+        TU1 = 10; thr[TU1] = 12; ch[8] = TU1; name[8]="TU1"
+        TU2 = 13; thr[TU2] = 12; ch[9] = TU2; name[9]="TU2"
+        TU3 = 8; thr[TU3] = 12; ch[10] = TU3; name[10]="TU3"
+        TU4 = 7; thr[TU4] = 12; ch[11] = TU4; name[11]="TU4"
+        TU5 = 4; thr[TU5] = 12; ch[12] = TU5; name[12]="TU5"
+        TU6 = 14; thr[2] = 12; ch[13] = 2; name[13] = "TU6"
+        TU7 = 1; thr[2] = 12; ch[14] = 2; name[14] = "TU7"
+        TU8 = 2; thr[2] =12; ch[15] = 2; name[15] = "TU8"
+        
+'''
+class Hist:
+        def __init__(self, dat_name = '2023-12-06-AHBGOC-Bi207-07.dat' ):
+                self.dat_name = dat_name
+                self.my_path2irena = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/data/'
+                self.chunksize = 2e6
+                # conversion to mV because, i of some adc conversion
+                self.mV = 14000
+                # name: channel
+                self.upper_trigger = {
+                'TO2': 3,
+                'TO5': 5,
+                'TO9': 6,
+                'TO13': 9,
+                'TO16': 11,
+                'TO11': 15,
+                'TO7': 0,
+                'TO12': 12
+                }
+
+                self.lower_trigger = {
+                'TU2': 10,
+                'TU5': 13,
+                'TU9': 8,
+                'TU13': 7,
+                'TU16': 4,
+                'TU11': 14,
+                'TU7': 1, 
+                'TU12': 2
+                }
+                self.calibration_factor_b1 = 0.652
+                self.calibration_factor_b2 = 0.668
+                
+                #self.create_EI_EI()
+                # Bezeichnung passt, alles abgeglichen mit meinen Notizen, die ich mir gemacht habe. Nochmal beim Umbau überprüfen und dann passt das.
+                self.bgo = {'B1':17,'B2':16}
+                #the channel corresponding to the name
+                self.minV = -100
+                self.maxV = 3500
+                self.resV = 0.838214/2
+                self.bins = np.arange(self.minV, self.maxV, self.resV)
+                #self.get_time()
+                self.create_hist_1D()
+
+        def to_mV(self, chunk):
+                # remove non EI lines   
+                for column in chunk.columns[5:]:
+                        if 'mV' in column:
+                                chunk[column] = chunk[column] / self.mV 
+                return chunk
+        
+        def init_results(self):
+                self.results = []
+                self.namelist = []
+                for i in self.upper_trigger:
+                        for j in self.lower_trigger:
+                                name = i+'_'+j
+                                # name, BGO1, BGO2, BGO1 enerrgy deposit, BGO2 energy deposit, triggerdioden counts, summed energy deposit
+                                self.namelist.append(name)
+                                self.results.append([name, np.zeros(len(self.bins)-1), np.zeros(len(self.bins)-1),0,0,0, np.zeros(len(self.bins)-1)])
+                                if i == 'TO8' and j == 'TU8':
+                                        break
+
+        def search_name_index_in_results(self,results, search_name):
+                for index, result in enumerate(results):
+                        if result[0] == search_name:
+                                return index
+                
+        def EI_columns(self):
+                 # create header for the data
+                self.header = ['Datenprodukt','Temperatur', 'timestamp', 'L2_trigger', 'timediff']
+                for i in range(0, 18):
+                        self.header.append('mV_'+str(i))
+                        self.header.append('trigger_'+str(i))
+                        self.header.append('phase_'+str(i))             
+                      
+        def create_hist_1D(self):
+                b1_hist = np.zeros(len(self.bins)-1)
+                b2_hist = np.zeros(len(self.bins)-1)
+                self.init_results()  
+                #print(np.shape(self.results))     
+                # which of the top diodes triggered, only EI lines
+                self.csv_reader = pd.read_csv(self.my_path2irena+self.dat_name[:-4]+'.EI.EI_data', chunksize=self.chunksize, delim_whitespace=True, header=None, engine ='c') 
+                for chunk in self.csv_reader:
+                        print('in chunk')
+                        self.EI_columns()
+                        chunk.columns = self.header
+                        chunk = self.to_mV(chunk)
+                        #chunk = self.B1_B2_triggered(chunk)
+                        # iterate over upper and lower trigger diodes
+                        # either 1 or 2 of the BGO's triggered 
+                        for i in self.upper_trigger:
+                                for j in self.lower_trigger:
+                                        # 2 diodes triggered
+                                        mask = (chunk['mV_'+str(self.upper_trigger[i])] > 12) & (chunk['mV_'+str(self.lower_trigger[j])] > 12) 
+                                        # bgo triggered as well
+                                        mask = mask & (chunk['mV_17'] < 12) & (chunk['mV_16'] < 12)
+                                        new_chunk = chunk[mask]
+                                        name = i+'_'+j 
+                                        hist_B1, bins_B1 = np.histogram(new_chunk['mV_17'], bins=self.bins)
+                                        hist_B2, bins_B2 = np.histogram(new_chunk['mV_16'], bins=self.bins)
+                                        #self.results[self.search_name_index_in_results(self.results, name)][1] += hist_B1
+                                        #self.results[self.search_name_index_in_results(self.results, name)][2] += hist_B2
+                                        #self.results[self.search_name_index_in_results(self.results, name)][3] += np.sum(new_chunk['mV_17'])
+                                        #self.results[self.search_name_index_in_results(self.results, name)][4] += np.sum(new_chunk['mV_16'])
+                                        self.results[self.search_name_index_in_results(self.results, name)][5] += len(new_chunk['mV_17'])
+                                        # hole energy, save in keV_Si, check whether it is the correct calibration factor
+                                        #hist_sum, bins_sum = np.histogram(new_chunk['mV_17']/self.calibration_factor_b1+new_chunk['mV_16']/self.calibration_factor_b2, bins=self.bins)
+                                        #self.results[self.search_name_index_in_results(self.results, name)][6] += hist_sum
+                # middle of the bins
+                self.bins = bins_B1[:-1] + (bins_B1[1] - bins_B1[0])/2
+                self.plot_missed_counts()
+                #self.write_file()
+        def plot_missed_counts(self):
+                print('start plotting')
+                # extract all labels from self.results
+                labels = []
+                counts = []
+                for i in self.results:
+                        labels.append(i[0])
+                        counts.append(i[5])
+                # sort labels and ratios alphabetically
+                labels, counts = zip(*sorted(zip(labels, counts)))
+                # put in pandasdataframe and save a csv file
+                df = pd.DataFrame({'Coincidence': labels, 'Counts': counts})
+                df.to_csv('missed_counts_12mV.csv', index=False)
+                # plot missed counts
+                plt.figure(figsize=(20, 10))
+                plt.bar(labels, counts)
+                plt.grid()
+                plt.title('Counts without BGO trigger but with 2 trigger diodes')
+                plt.ylabel('Counts')
+                plt.xticks(rotation=90)
+                plt.savefig('missed_counts_12mV.pdf', bbox_inches='tight')
+
+                plt.show()
+
+        def write_file(self):
+                print('write file')
+                # back to working directory
+                os.chdir('/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/')
+                # Ensure the 'Results' folder exists
+                main_folder = f'Results{self.dat_name[:-4]}'
+                muon_folder = os.path.join(main_folder,'1d_hist_muon_check_dark_counts')
+                os.makedirs(muon_folder, exist_ok=True)
+                # Create a file for each coincidence
+                for i in self.results:
+                        name = i[0]
+                        B1 = i[1]
+                        B2 = i[2]
+                        counts = i[5]
+                        sum_b1_b2 = i[6]
+                        file_path = os.path.join(muon_folder, f'{name}.txt')
+                        # List of comments
+                        comments = [
+                        f'# measurement time in s: {self.measurement_time[:-1]}',
+                        f'# BGO1 energy deposit in mV: {i[3]}',
+                        f'# BGO2 energy deposit in mV: {i[4]}',
+                        f'# counts: {counts}\n'
+                        ]
+
+                        # Save data to CSV file with header comments
+                        header_comments = '\n'.join(comments)
+                        np.savetxt(file_path, np.column_stack((self.bins, B1, B2, sum_b1_b2)), delimiter=',',
+                                header=f'#BGO1,BGO2', comments=header_comments)
+
+        def create_EI_EI(self):
+                os.chdir(self.my_path2irena)
+                os.system(f'rm {self.my_path2irena+self.dat_name[:-4]}.EI.EI_data')
+                os.system(f'awk -f filter.awk {self.my_path2irena+self.dat_name[:-4]}.EI')
+        
+        def get_time(self):
+                os.chdir(self.my_path2irena)
+                self.measurement_time = os.popen(f'awk -f Itime.awk {self.my_path2irena+self.dat_name[:-4]}.EI').read()
+
+
+source_file = '2023-12-21-AHBGOC-14-langzeit.dat'
+Hist(dat_name=source_file)
+
+
+
+        

create_1d_hist_muon.py

@@ -143,7 +143,7 @@ class Hist:
                                         # 2 diodes triggered
                                         mask = (chunk['mV_'+str(self.upper_trigger[i])] > 12) & (chunk['mV_'+str(self.lower_trigger[j])] > 12) 
                                         # bgo triggered as well
-                                        mask = mask & (chunk['mV_17'] > 12) & (chunk['mV_16'] > 12)
+                                        #mask = mask & (chunk['mV_17'] > 12) & (chunk['mV_16'] > 12)
                                         new_chunk = chunk[mask]
                                         name = i+'_'+j +'m'
                                         hist_B1, bins_B1 = np.histogram(new_chunk['mV_17'], bins=self.bins)
@@ -178,7 +178,7 @@ class Hist:
                         file_path = os.path.join(muon_folder, f'{name}.txt')
                         # List of comments
                         comments = [
-                        f'# measurement time in s: {self.measurement_time[:-2]}',
+                        f'# measurement time in s: {self.measurement_time[:-1]}',
                         f'# BGO1 energy deposit in mV: {i[3]}',
                         f'# BGO2 energy deposit in mV: {i[4]}',
                         f'# counts: {counts}\n'

create_2d_hist_muon.py

@@ -38,7 +38,7 @@ class Hist:
                 self.dat_name = '2023-12-21-AHBGOC-14-langzeit.EI.EI_data'
                 self.saving_path = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit'
                 self.my_path2irena = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/data/'
-                self.chunksize = 3e6
+                self.chunksize = 2e6
                 self.calibration_factor_b1 = 0.652
                 self.calibration_factor_b2 = 0.668
                 # to mV
@@ -65,6 +65,7 @@ class Hist:
                 'TU7': 1,
                 'TU12': 2
                 }
+                # this is still the old BGO mapping
                 self.bgo = {'B1':17,'B2':16}
                 self.minV = -100
                 self.maxV = 3500

create_muon_hist.py

@@ -0,0 +1,197 @@
+import numpy as np
+import pandas as pd
+import os
+import matplotlib.pyplot as plt
+import tqdm
+plt.rcParams.update({'font.size': 25})
+
+'''
+This script outputs the xray histogramms
+
+
+Kanalzuordnung:
+        TO1 = 0; thr[TO1] = 12; ch[0] = TO1; name[0]="TO1"
+        TU1 = 1; thr[TU1] = 12; ch[1] = TU1; name[1]="TU1"
+        TU2 = 2; thr[TU2] = 12; ch[2] = TU2; name[2]="TU2"
+        TO2 = 3; thr[TO2] = 12; ch[3] = TO2; name[3]="TO2"
+        TU3 = 4; thr[TU3] = 12; ch[4] = TU3; name[4]="TU3"
+        TO3 = 5; thr[TO3] = 12; ch[5] = TO3; name[5]="TO3"
+        TO4 = 6; thr[TO4] = 12; ch[6] = TO4; name[6]="TO4"
+        TU4 = 7; thr[TU4] = 12; ch[7] = TU4; name[7]="TU4"
+        TU5 = 8; thr[TU5] = 12; ch[8] = TU5; name[8]="TU5"
+        TO5 = 9; thr[TO5] = 12; ch[9] = TO5; name[9]="TO5"
+        TU7 = 10; thr[TU7] = 12; ch[10] = TU7; name[10] = "TU7"
+        TO6 = 11; thr[11] = 12; ch[11] = TO6; name[11] = "TO6"
+        TO7 = 12; thr[12] = 12; ch[12] = TO7; name[12] = "TO7"
+        TU6 = 13; thr[13] = 12; ch[13] = 13; name[13] = "TU6"
+        TU8 = 14; thr[14] = 12; ch[14] = 14; name[14] = "TU8"
+        TO8 = 15; thr[15] =12; ch[15] = 15; name[15] = "TO8"
+        B2 = 16; thr[16] =12; ch[16] = 16; name[16] = "B2"
+        B1 = 17; thr[17] = 12; ch[17] = 17; name[17] = "B1" 
+
+neue Kanalzuordnung:
+        oben:
+        TO1 = 3; thr[TO1] = 12; ch[0] = TO1; name[0]="TO1" 
+        TO2 = 5; thr[TO2] = 12; ch[1] = TO2; name[1]="TO2"
+        TO3 = 6; thr[TO3] = 12; ch[2] = TO3; name[2]="TO3"
+        TO4 = 9; thr[TO4] = 12; ch[3] = TO4; name[3]="TO4"
+        TO5 = 11; thr[TO5] = 12; ch[4] = TO5; name[4]="TO5"
+        TO6 = 15; thr[15] = 12; ch[5] = TO6; name[5] = "TO6"
+        TO7 = 0; thr[0] = 12; ch[6] = TO7; name[6] = "TO7"
+        TO8 = 12; thr[12] =12; ch[7] = 12; name[7] = "TO8"
+        unten:
+        TU1 = 10; thr[TU1] = 12; ch[8] = TU1; name[8]="TU1"
+        TU2 = 13; thr[TU2] = 12; ch[9] = TU2; name[9]="TU2"
+        TU3 = 8; thr[TU3] = 12; ch[10] = TU3; name[10]="TU3"
+        TU4 = 7; thr[TU4] = 12; ch[11] = TU4; name[11]="TU4"
+        TU5 = 4; thr[TU5] = 12; ch[12] = TU5; name[12]="TU5"
+        TU6 = 14; thr[2] = 12; ch[13] = 2; name[13] = "TU6"
+        TU7 = 1; thr[2] = 12; ch[14] = 2; name[14] = "TU7"
+        TU8 = 2; thr[2] =12; ch[15] = 2; name[15] = "TU8"
+        
+'''
+class Hist:
+        def __init__(self, dat_name = '2023-12-21-AHBGOC-14-langzeit.dat' ):
+                self.rem_path = '/data/asterix/athena/arm/irena/ahbgo/'
+                self.dat_name = dat_name
+                self.my_path2irena = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/data/'
+                self.chunksize = 2e6
+                # to mV
+                self.mV = 14000
+                # name: channel
+                self.upper_trigger = {
+                'TO1': 3,
+                'TO2': 5,
+                'TO3': 6,
+                'TO4': 9,
+                'TO5': 11,
+                'TO6': 15,
+                'TO7': 0,
+                'TO8': 12
+                }
+
+                self.lower_trigger = {
+                'TU1': 10,
+                'TU2': 13,
+                'TU3': 8,
+                'TU4': 7,
+                'TU5': 4,
+                'TU6': 14,
+                'TU7': 1,
+                'TU8': 2
+                }
+
+                #self.create_EI_EI()
+                self.bgo = {'B1':17,'B2':16}
+                #the channel corresponding to the name
+                self.minV = -100
+                self.maxV = 3500
+                self.resV = 0.838214/2
+                self.bins = np.arange(self.minV, self.maxV, self.resV)
+                print(self.bins)
+                self.get_time()
+                self.create_hist_1D()
+
+        def to_mV(self, chunk):
+                # remove non EI lines   
+                for column in chunk.columns[5:]:
+                        if 'mV' in column:
+                                chunk[column] = chunk[column] / self.mV 
+                return chunk
+        
+        def init_results(self):
+                self.results = []
+                self.namelist = []
+                for i in {**self.upper_trigger, **self.lower_trigger, **self.bgo}:
+                        # name, BGO1, BGO2, BGO1 enerrgy deposit, BGO2 energy deposit
+                        self.namelist.append(i+'_xray')
+                        self.results.append([i+'_xray', np.zeros(len(self.bins)-1)])
+
+
+        def search_name_index_in_results(self,results, search_name):
+                for index, result in enumerate(results):
+                        if result[0] == search_name:
+                                return index
+                
+        def EI_columns(self):
+                 # create header for the data
+                self.header = ['Datenprodukt','Temperatur', 'timestamp', 'L2_trigger', 'timediff']
+                for i in range(0, 18):
+                        self.header.append('mV_'+str(i))
+                        self.header.append('trigger_'+str(i))
+                        self.header.append('phase_'+str(i))             
+                      
+        def create_hist_1D(self):
+                print('create hist')
+                # which of the top diodes triggered, only EI lines
+                print('path',self.my_path2irena+self.dat_name[:-4]+'.EI.EI_data')
+                self.csv_reader = pd.read_csv(self.my_path2irena+self.dat_name[:-4]+'.EI.EI_data', chunksize=self.chunksize, delim_whitespace=True, header=None, engine='c') 
+                self.init_results()       
+                for chunk in self.csv_reader:
+                        self.EI_columns()
+                        chunk.columns = self.header
+                        chunk = self.to_mV(chunk)
+                        combined_trigger = {**self.upper_trigger, **self.lower_trigger, **self.bgo}
+                        for i in combined_trigger:
+                            if i != 'B1' and i != 'B2':
+                                mask = (chunk['mV_'+str(combined_trigger[i])] > 20)
+                                mask = mask  & ((chunk['mV_17'] > 20) & (chunk['mV_16'] > 20))
+                                new_chunk = chunk[mask]
+                                name = i+'_xray'
+                                hist, bins = np.histogram(new_chunk['mV_'+str(combined_trigger[i])], bins=self.bins)
+                                self.results[self.search_name_index_in_results(self.results, name)][1] += hist
+                            else:
+                                # for relevant ran
+                                mask = ((chunk['mV_17'] > 20) & (chunk['mV_16'] > 20))
+                                new_chunk = chunk[mask]
+                                name = i+'_xray'
+                                hist, bins = np.histogram(new_chunk['mV_'+str(combined_trigger[i])], bins=self.bins)
+                                self.results[self.search_name_index_in_results(self.results, name)][1] += hist
+                # save middle of bins
+                self.bins = bins[:-1] + self.resV/2
+                self.write_file()
+
+        def write_file(self):
+                print('write file')
+                # back to working directory
+                os.chdir('/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/')
+                # Ensure the 'Results' folder exists
+                main_folder = f'Results{self.dat_name[:-4]}'
+                xray_folder = os.path.join(main_folder,'1d_hist_muon_no_coinc')
+                os.makedirs(xray_folder, exist_ok=True)
+                # Create a file for each coincidence
+                for i in self.results:
+                        name = i[0]
+                        counts = i[1]
+                        file_path = os.path.join(xray_folder, f'{name}.txt')
+                        # List of comments
+                        comments = [
+                        f'# measurement time in s: {self.measurement_time}']
+
+                        # Save data to CSV file with header comments
+                        np.savetxt(file_path, np.column_stack((self.bins, counts)), delimiter=',',
+                                header=f'BGO1,BGO2', comments='\n'.join(comments))
+
+
+        def create_EI_EI(self):
+                print('create EI_EI')
+                os.chdir(self.my_path2irena)
+                os.system(f'rm {self.my_path2irena+self.dat_name[:-4]}.EI.EI_data')
+                print('create EI_EI_start')
+                os.system(f'awk -f filter.awk {self.my_path2irena+self.dat_name[:-4]}.EI')
+                print('create EI_EI_done')
+        
+        def get_time(self):
+                print('get time')
+                os.chdir(self.my_path2irena)
+                self.measurement_time = os.popen(f'awk -f Itime.awk {self.my_path2irena+self.dat_name[:-4]}.EI').read()
+        
+
+
+
+source_file = '2023-12-21-AHBGOC-14-langzeit.dat'
+Hist(dat_name=source_file)
+
+
+
+        

create_xray_hist.py

@@ -51,7 +51,7 @@ neue Kanalzuordnung:
         
 '''
 class Hist:
-        def __init__(self, dat_name = '2023-12-21-AHBGOC-14-langzeit.dat' ):
+        def __init__(self, dat_name = '2023-12-11-AHBGOC-Bi207-08.dat' ):
                 self.rem_path = '/data/asterix/athena/arm/irena/ahbgo/'
                 self.dat_name = dat_name
                 self.my_path2irena = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/data/'
@@ -82,6 +82,7 @@ class Hist:
                 }
 
                 #self.create_EI_EI()
+                # this is the old mapping
                 self.bgo = {'B1':17,'B2':16}
                 #the channel corresponding to the name
                 self.minV = -100
@@ -140,9 +141,16 @@ class Hist:
                                 name = i+'_xray'
                                 hist, bins = np.histogram(new_chunk['mV_'+str(combined_trigger[i])], bins=self.bins)
                                 self.results[self.search_name_index_in_results(self.results, name)][1] += hist
-                            else:
+                            if i == 'B1':
                                 # for relevant ran
-                                mask = ((chunk['mV_17'] > 12) ^ (chunk['mV_16'] > 12))
+                                mask = ((chunk['mV_17'] > 12) & (chunk['mV_16'] < 12))
+                                new_chunk = chunk[mask]
+                                name = i+'_xray'
+                                hist, bins = np.histogram(new_chunk['mV_'+str(combined_trigger[i])], bins=self.bins)
+                                self.results[self.search_name_index_in_results(self.results, name)][1] += hist
+                            if i == 'B2':
+                                # for relevant ran
+                                mask = ((chunk['mV_17'] < 12) & (chunk['mV_16'] > 12))
                                 new_chunk = chunk[mask]
                                 name = i+'_xray'
                                 hist, bins = np.histogram(new_chunk['mV_'+str(combined_trigger[i])], bins=self.bins)
@@ -190,7 +198,7 @@ class Hist:
 
 
 source_file = '2023-12-21-AHBGOC-14-langzeit.dat'
-Hist(dat_name=source_file)
+Hist()
 
 
 

cut_aftertime_EI.py

@@ -0,0 +1,15 @@
+from datetime import datetime
+import pandas as pd
+
+
+class CutAfterTime:
+    def __init__(self) -> None:
+        self.path = ''
+        self.new_EI_data = []
+    
+    def iterate_over(self):
+        # read chunkwise with pandas
+        df_contructor = pd.read_csv(self.path, chunksize=1e5, delimiter=',')
+        for df in df_contructor:
+            pass
+        

getamountofcounts.py

@@ -0,0 +1,43 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import os
+plt.rcParams.update({'font.size': 18})
+
+class nuumberofcoinc:
+    def __init__(self) -> None:
+        self.path = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/1d_hist_muon/'
+        self.names = []
+        self.counts = []    
+        self.get_names()
+        self.get_counts()
+        self.plot_bar_chart()
+    
+    def get_names(self):
+        self.file_names = os.listdir(self.path)
+    
+    def get_counts(self):
+        for i in range(len(self.file_names)):
+            data = np.loadtxt(self.path + self.file_names[i], comments='#', delimiter=',')
+            self.names.append(self.file_names[i][:-5])
+            self.counts.append(np.sum(data[:,1]))
+        # create pandas dataframe and save as csv
+        df = pd.DataFrame({'Coincidence': self.names, 'Counts': self.counts})
+        df.to_csv('counts12mV.csv', index=False)
+
+    
+    def plot_bar_chart(self):
+        # sort the names and counts
+        labels, counts = zip(*sorted(zip(self.names, self.counts)))
+        # plot
+        plt.figure(figsize=(20,10))
+        plt.bar(labels, counts)
+        plt.xticks(rotation=90)
+        plt.xlabel('Coincidence')
+        plt.ylabel('Counts')
+        plt.title('Number of coincidences with 2 trigger diodes')
+        plt.grid()
+        plt.savefig('counts12mV.pdf', bbox_inches='tight')
+        plt.show()
+
+nuumberofcoinc()

investigate_yield.py

@@ -0,0 +1,39 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import os
+plt.rcParams.update({'font.size': 18})
+class CheckYield:
+    def __init__(self):
+        self.file_path_only_trigger = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/1d_hist_muon_only_trigger_diodes/'
+        self.file_path_all_trigger = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/1d_hist_muon/'
+        self.file_names = os.listdir(self.file_path_only_trigger)
+        self.chech_ratio_counts_and_plot()
+
+    def chech_ratio_counts_and_plot(self):
+        ratios = []
+        labels = []
+        plt.figure(figsize=(20, 10))
+        for i in self.file_names:
+            labels.append(i[:-5])
+            data_only_trigger = np.loadtxt(self.file_path_only_trigger + i, comments='#', delimiter=',')
+            data_all_trigger = np.loadtxt(self.file_path_all_trigger + i, comments='#', delimiter=',')
+            hist_only_trigger = data_only_trigger[:,1]
+            hist_all_trigger = data_all_trigger[:,1]
+            ratios.append(abs(np.sum(hist_only_trigger)-np.sum(hist_all_trigger)))
+
+            # bar plot of ratios
+        # how much percent higher
+        print('Overall percent higher: ', np.mean(ratios))
+        # sort labels and ratios alphabetically
+        labels, ratios = zip(*sorted(zip(labels, ratios)))
+        plt.bar(labels, ratios)
+        plt.xticks(rotation=90)
+        plt.title('Countdifference 1d histogram with only trigger diodes and trigger diodes plus BGO diodes')
+        plt.ylabel('Countdiffference')
+        #plt.yscale('log')
+        plt.grid()
+        plt.savefig('difference_counts.pdf', bbox_inches='tight')
+        plt.show()
+
+CheckYield()

orthogonal_mapping_1d.py

@@ -6,11 +6,14 @@ from scipy.optimize import curve_fit
 from matplotlib.ticker import FixedLocator
 import matplotlib.image as image
 from matplotlib.offsetbox import (OffsetImage, AnnotationBbox)
+from round_DIN1333 import  round_DIN1333
+from decimal import Decimal, getcontext
 import math
 class plot_1d_hist:
     def __init__(self) -> None:
-        self.path = 'Results2023-12-21-AHBGOC-14-langzeit/1d_hist_muon'
+        self.path = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/1d_hist_muon'
         self.file_names = self.file_names()
+        self.silizium = 3.6e-3
         self.mapping = pd.DataFrame({
             # vertikale linie von diode zu diode
             2: [0,40], 
@@ -39,6 +42,7 @@ class plot_1d_hist:
         self.energy_photon_keV = 3.6e-3
         self.used = [2,5,9,13,16,11,7,12]
         self.right = 7
+        self.resV = 0.838214/2
         self.create_plot()
         self.mpvs_b1 = []
         self.mpvs_b2 = []
@@ -71,52 +75,65 @@ class plot_1d_hist:
             # load in if 2 times the same number in file 
             file = [os.path.basename(file_path) for file_path in self.file_names if os.path.basename(file_path).startswith(f'TO{i}_TU{i}')]
             # load file
+            print('path', self.path+'/'+file[0])
             data = np.loadtxt(self.path+'/'+file[0], delimiter=',', comments='#')
             # lese Messzeit aus erste Zeile ein
             with open(self.path+'/'+file[0], 'r') as file:
                 measurement_time =  int(next(file).split(':')[-1].strip()) / 60 / 60
-            B2 = data[:, 1] / measurement_time
-            B1 = data[:, 2] / measurement_time
-            B1_B2 = data[:, 3] / measurement_time
+                self.measurement_time = measurement_time
+            B2 = data[:, 1] / measurement_time / self.resV
+            B1 = data[:, 2] / measurement_time / self.resV
+            B1_B2 = data[:, 3] / measurement_time / self.resV
 
             bins_b1 = data[:, 0] / self.calibration_factor_b2
             bins_b2 = data[:, 0] / self.calibration_factor_b1
-            bins_b1_b2 = data[:, 0] / self.calibration_factor_b1
+            bins_b1_b2 = data[:, 0]
             # fit landau
             x_1, y_1, mpv_1, e_error_1 = self.fit_landau(bins_b1,B1)
             x_2, y_2, mpv_2, e_error_2 = self.fit_landau(bins_b2,B2)
             x_b1b2, y_b1b2, mpv_b1b2, e_error_b1b2 = self.fit_landau(bins_b1_b2,B1_B2, sum = True)
 
+            mpv_1_round, e_error_1_round = round_DIN1333(mpv_1/self.silizium/100000, e_error_1/self.silizium/100000)
+            mpv_2_round, e_error_2_round = round_DIN1333(mpv_2/self.silizium/100000, e_error_2/self.silizium/100000)
+            mpv_b1b2_round, e_error_b1b2_round = round_DIN1333(mpv_b1b2/self.silizium/100000, e_error_b1b2/self.silizium/100000)
+            mpv_1_round = self.multiply_with_error(mpv_1_round)
+            mpv_2_round = self.multiply_with_error(mpv_2_round)
+            mpv_b1b2_round = self.multiply_with_error(mpv_b1b2_round)
+            e_error_1_round = self.multiply_with_error(e_error_1_round)
+            e_error_2_round = self.multiply_with_error(e_error_2_round)
+            e_error_b1b2_round = self.multiply_with_error(e_error_b1b2_round)
+
 
             # pass 
             name = fr'$m_{{\text{{TO{i}TU{i}B1}}}}$ = '
-            ax.step(bins_b1, B1, label=f'B1 mpv: m = '+ f'({round(mpv_1,2)}$\pm${math.ceil(e_error_1*100)/100})'+ ' $keV_{Si}$', where='mid')
+
+            ax.step(bins_b1/self.silizium, B1, label=f'B1 mpv: m = '+ f'{mpv_1_round}$\pm${e_error_1_round}', where='mid')
             name = f'TO{i}_TU{i}_B2'
-            ax.step(bins_b2, B2, label='B2 mpv: m = ' + f'({round(mpv_2,2)}$\pm${math.ceil(e_error_2*100)/100})'+ ' $keV_{Si}$', where='mid')
+            ax.step(bins_b2/self.silizium, B2, label='B2 mpv: m = ' + f'{mpv_2_round}$\pm${e_error_2_round}', where='mid')
             name = f'TO{i}_TU{i}_B1+B2'
-            ax.step(bins_b1_b2, B1_B2, label='B1+B2 mpv: m = ' + f'({round(mpv_b1b2,2)}$\pm${math.ceil(e_error_b1b2*100)/100})'+ ' $keV_{Si}$', where='mid')
-            plt.plot(x_1, y_1, linestyle='--', color='black')
-            plt.plot(x_2, y_2, linestyle='--', color='black')
-            plt.plot(x_b1b2, y_b1b2, linestyle='--', color='black')
+            ax.step(bins_b1_b2/self.silizium, B1_B2, label='B1+B2 mpv: m = ' + f'{mpv_b1b2_round}$\pm${e_error_b1b2_round}', where='mid')
+            plt.plot(x_1/self.silizium, y_1, linestyle='--', color='black')
+            plt.plot(x_2/self.silizium, y_2, linestyle='--', color='black')
+            plt.plot(x_b1b2/self.silizium, y_b1b2, linestyle='--', color='black')
     
             # Deactivate x-axis labels
             if i  == self.right:
-                ax.set_xticks([100,200,300])
+                ax.set_xticks([25000,50000,75000])
                 ax.grid(True)
             else: 
                 # Adjust x-axis ticks for compactness
-                ax.set_xticks([100,200,300,400])
+                ax.set_xticks([25000,50000,75000, 100000])
                 ax.grid(True)
             
             ax.set_yscale('log')
             # y.axis label
             x_ticks = [0.001,0.1]
             ax.yaxis.set_major_locator(FixedLocator(x_ticks))
-            ax.set_ylim(bottom = 0.9/ measurement_time)
-            ax.set_xlim(0, 500)
+            ax.set_ylim(bottom = 0.9/ measurement_time, top = 1)
+            ax.set_xlim(0, 500/self.silizium)
             ax.legend(loc='upper right', bbox_to_anchor=(1.2, 1.3))
         # Adding a text block
-        text_block = 'x-axis: Energydeposit diode per diode in $keV_{Si}$\n y-axis: Number of events in 1/h/mV'
+        text_block = 'x-axis: emitted photoelectrons\n y-axis: number of events in 1/h/mV'
         bbox_props = dict(boxstyle='round', facecolor='white', edgecolor='black', linewidth=1)
         plt.text(0.25, 0.25, text_block, transform=fig.transFigure, fontsize=12, va='center', ha='center', bbox=bbox_props)
 
@@ -133,8 +150,8 @@ class plot_1d_hist:
         plt.show()
 
     def fit_landau(self, x_data, y_data, sum = False):
-        mid_cut_low = 20
-        mid_cut_high = 200
+        mid_cut_low = 30
+        mid_cut_high = 130
         if sum:
             mid_cut_low = 80
             mid_cut_high = 400
@@ -143,12 +160,13 @@ class plot_1d_hist:
         y_s = []
         mpv_s = []
         e_errors = []
-        for lower_cut in range(0,10,2):
+        for lower_cut in range(0,20,2):
             for higher_cut in range(0,50,10):
                 try:   
                     mask = (x_data > mid_cut_low+lower_cut) & (x_data < mid_cut_high+higher_cut)
                     new_x, y, mpv, e_error = self.fit_function_landau(x_data[mask], y_data[mask])
                     error = np.sum(np.abs(y_data[mask][1:] - y))
+                    print(error)
                     x_s.append(new_x)
                     y_s.append(y)
                     mpv_s.append(mpv)
@@ -162,14 +180,13 @@ class plot_1d_hist:
 
     def fit_function_landau(self, x_data, y_data):
         # fit on middle of the bins and not on the edges
-        x_data = x_data[:-1] + (x_data[1] - x_data[0]) / 2
+        x_data = x_data[:-1]
         y_data = y_data[1:]        
         # First iteration
-        sigma = np.sqrt(y_data+1)
-        p0 = [np.max(y_data), x_data[np.argmax(y_data)]-10, 10, np.max(y_data), 100]
-        popt, pcov = curve_fit(self.lan, x_data, y_data, p0=p0, maxfev = 1000, sigma=sigma, absolute_sigma=True)
-
-        a_error, e_error, s_error, ab_error, sb_error = np.sqrt(np.diag(pcov))
+        sigma = np.sqrt(y_data+1/self.resV/self.measurement_time)
+        p0 = [np.max(y_data), x_data[np.argmax(y_data)], 10]
+        popt, pcov = curve_fit(self.lan, x_data, y_data, p0=p0, maxfev = 2000, sigma=sigma, absolute_sigma=True)
+        a_error, e_error, s_error = np.sqrt(np.diag(pcov))
         #new_x = np.linspace(x_data[0], x_data[-1], 100)
         y = self.lan(x_data, *popt)
         mpv = popt[1]
@@ -187,7 +204,7 @@ class plot_1d_hist:
         return a*np.exp((e-x)/s)
 
     
-    def lan(self,x, a,e,s,ab,sb):
+    def lan(self,x, a,e,s):
         '''
         a = amplitude
         e = position
@@ -195,6 +212,15 @@ class plot_1d_hist:
         ab = amplitude bg
         sb = width bg
         '''
-        return self.mips(x, a,e,s) + self.bg(x, ab, e, sb)
+        return self.mips(x, a,e,s)
+    def multiply_with_error(self,value):
+        # Set precision to a sufficiently high value
+        getcontext().prec = 50  # You can adjust the precision as needed
+        # Convert values to Decimal
+        value_decimal = Decimal(str(value))
+        factor_decimal = Decimal(str(100000))
+        # Perform multiplication with arbitrary precision
+        new_value = value_decimal * factor_decimal
+        return float(str(new_value).rstrip('0').rstrip('.'))
 
 plot_1d_hist()

orthogonal_mapping_2d.py

@@ -8,10 +8,11 @@ import matplotlib.image as image
 from matplotlib.offsetbox import (OffsetImage, AnnotationBbox)
 import math
 from matplotlib import cm
-
+from decimal import Decimal, getcontext
+from round_DIN1333 import  round_DIN1333
 class plot_1d_hist:
     def __init__(self) -> None:
-        self.path = 'Results2023-12-21-AHBGOC-14-langzeit/2dhistmuons'
+        self.path = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/2dhistmuons'
         self.file_names = self.file_names()
         self.mapping = pd.DataFrame({
             # vertikale linie von diode zu diode
@@ -38,7 +39,7 @@ class plot_1d_hist:
         })
         self.calibration_factor_b1 = 0.652
         self.calibration_factor_b2 = 0.668
-        self.energy_photon_keV = 3.6e-3
+        self.silizium = 3.6e-3
         self.used = [2,5,9,13,16,11,7,12]
         self.right = 7
         self.minV = -100
@@ -49,6 +50,15 @@ class plot_1d_hist:
         self.bins = np.arange(self.minV, self.maxV, self.resV)
         self.create_plot()
 
+    def multiply_with_error(self,value, multiplier):
+        # Set precision to a sufficiently high value
+        getcontext().prec = 50  # You can adjust the precision as needed
+        # Convert values to Decimal
+        value_decimal = Decimal(str(value))
+        factor_decimal = Decimal(str(multiplier))
+        # Perform multiplication with arbitrary precision
+        new_value = value_decimal * factor_decimal
+        return float(str(new_value).rstrip('0').rstrip('.'))
 
     def file_names(self):
         file_names = []
@@ -82,39 +92,38 @@ class plot_1d_hist:
             self.data = np.loadtxt(self.path+'/'+file[0], delimiter=',', comments='#')
             # lese Messzeit aus erste Zeile ein
             self.bins_new = self.bins[(self.bins >= xlim[0] ) & (self.bins <= xlim[1]) & (self.bins >= ylim[0]) & (self.bins <= ylim[1])]
-            self.bins_new_b1 = self.bins_new / self.calibration_factor_b1
-            self.bins_new_b2 = self.bins_new / self.calibration_factor_b2
+            self.bins_new_b1 = self.bins_new / self.calibration_factor_b2
+            self.bins_new_b2 = self.bins_new / self.calibration_factor_b1
                         
             self.data = self.data[(self.bins[:-1] >= xlim[0]) & (self.bins[:-1] <= xlim[1]), :]
             self.data = self.data[:, (self.bins[:-1] >= ylim[0]) & (self.bins[:-1] <= ylim[1])]
-            im = ax.pcolormesh(self.bins_new_b1, self.bins_new_b2, self.data, cmap=cmap, shading='auto', vmin=0.6)
+            im = ax.pcolormesh(self.bins_new_b2/ self.silizium, self.bins_new_b1/ self.silizium, self.data, cmap=cmap, shading='auto', vmin=0.6)
             x = y = np.arange(xlim[0],xlim[1], 1)
-            ax.plot(x, y, color='black',linewidth=2, linestyle='--',label='line with slope 1 and intercept 0')
+            ax.plot(x/self.silizium, y/self.silizium, color='black',linewidth=2, linestyle='--',label='line with slope 1 and intercept 0')
             self.line_data_from_matrice()
-            popt, pcov = curve_fit(self.line, self.x, self.y,p0=[1,0], method='lm')
+            popt, pcov = curve_fit(self.line, self.x, self.y,p0=[1,0])
             a_error = np.sqrt(pcov[0][0])
             b_error = np.sqrt(pcov[1][1])
-            ax.plot(x, self.line(x, *popt), color='red', linestyle='--', linewidth=2, label=f'line fit:\nslope s = ({round(popt[0],3)} ± {round(a_error,3)})\nintercept b = ({round(popt[1],1)} ± {round(b_error,1)})'+' $keV_{Si}$')
+            round_a_value, round_a_error = round_DIN1333(popt[0]/10, a_error/10) 
+            round_b_value, round_b_error = round_DIN1333(popt[1]/self.silizium/10000, b_error/self.silizium/10000)
+            round_a_value = self.multiply_with_error(round_a_value,10)
+            round_b_value = self.multiply_with_error(round_b_value,10000)
+            round_a_error = self.multiply_with_error(round_a_error,10)
+            round_b_error = self.multiply_with_error(round_b_error,10000)   
+            print('round_a_value: ', round_a_value, 'round_a_error: ', round_a_error, 'round_b_value: ', round_b_value, 'round_b_error: ', round_b_error)
+            ax.plot(x/ self.silizium, self.line(x, *popt)/ self.silizium, color='red', linestyle='--', linewidth=2, label=f'line fit:\nslope s = {round_a_value} ± {round_a_error}\nintercept b = {round_b_value} ± {round_b_error}')
 
             plt.grid()
             plt.legend(loc = 'upper left')
     
-
-            # Deactivate x-axis labels
-            if i  == self.right:
-                ax.set_xticks([75,100,125])
-                ax.grid(True)
-            else: 
-                # Adjust x-axis ticks for compactness
-                ax.set_xticks([75,100,125])
-                ax.grid(True)
-            ax.set_ylim(50,150)
-            ax.set_xlim(50,150)
+            ax.set_xticks([20000,25000,30000,35000])
+            ax.set_ylim(50/self.silizium,150/self.silizium)
+            ax.set_xlim(50/self.silizium,150/self.silizium)
 
             ax.legend(loc='upper right', bbox_to_anchor=(1.1, 1.3))
 
         # Adding a text block
-        text_block = 'x-axis: Energydeposit B2 in $keV_{Si}$\n y-axis: Energydeposit B1 in $keV_{Si}$'
+        text_block = 'x-axis: B2 emitted photoelectrons\ny-axis: B1 emitted photoelectrons'
         bbox_props = dict(boxstyle='round', facecolor='white', edgecolor='black', linewidth=1)
         plt.text(0.25, 0.25, text_block, transform=fig.transFigure, fontsize=12, va='center', ha='center', bbox=bbox_props)
 
@@ -127,8 +136,7 @@ class plot_1d_hist:
         fig.add_artist(ab)
                
         plt.suptitle('Overview of vertical coincidence measurements in top view', x=0.525, y=0.95, fontsize=20)
-        plt.savefig('Overview_vertical_coinc_top_view_2d_hist.png', dpi = 300)
-        plt.show()
+        plt.savefig('Overview_vertical_coinc_top_view_2d_hist.png', dpi = 400)
         
     def line_data_from_matrice(self):
         # i want to create histogram x and y data from the 2d histogram matrix
@@ -143,8 +151,8 @@ class plot_1d_hist:
                 for index_column, column in enumerate(row):
                         for times in range(int(self.data[index_row][index_column])):
                                 # x is B1 so 
-                                x.append(self.bins_new_b1[index_column])
-                                y.append(self.bins_new_b2[index_row])
+                                x.append(self.bins_new_b2[index_column])
+                                y.append(self.bins_new_b1[index_row])
 
         self.x = np.array(x)
         self.y = np.array(y)

plot_1d_hist_muon_no_coinc.py

@@ -0,0 +1,39 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+
+
+class Plot1DHist:
+    def __init__(self) -> None:
+        self.directory = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/1d_hist_muon_no_coinc'
+        self.file_names()   
+        self.plot_hist()
+    
+    def plot_hist(self):
+        counts_below = []
+        counts_top = []
+        for i in self.files:
+            data = np.loadtxt(i, delimiter=',', skiprows=2)
+            # mit kernel glätten
+            data[:,1] = np.convolve(data[:,1], np.ones(3), 'same') / 3
+            plt.step(data[:,0], data[:,1], where='mid', label=os.path.basename(i)[:-4])
+            # print counts
+            print(f'{os.path.basename(i)[:-4]}: {np.sum(data[:,1])}')
+            if 'U' in os.path.basename(i)[:-4]:
+                counts_below.append(np.sum(data[:,1]))
+            if 'O' in os.path.basename(i)[:-4]:
+                counts_top.append(np.sum(data[:,1]))
+        plt.legend()
+        plt.yscale('log')
+        plt.title(f'Muons of run 2023-12-21-AHBGOC-14-langzeit')
+        plt.show()
+        print('mean counts below: ', np.mean(counts_below))
+        print('mean counts top: ', np.mean(counts_top))
+        print('mean ratio: ', np.mean(counts_top)/np.mean(counts_below))
+    def file_names(self):
+        self.files = []
+        for file in os.listdir(self.directory):
+            if file.endswith(".txt"):
+                self.files.append(os.path.join(self.directory, file))
+
+Plot1DHist()

plot_1d_histogramm_muon.py

@@ -8,12 +8,15 @@ from scipy.optimize import curve_fit
 from matplotlib.ticker import MultipleLocator
 import re
 import math
+from decimal import Decimal, getcontext
 import matplotlib.image as image
 from matplotlib.offsetbox import (OffsetImage, AnnotationBbox)
 #ignore warnings
+
+from round_DIN1333 import round_DIN1333
 import warnings
 warnings.filterwarnings("ignore")
-plt.rcParams.update({'font.size': 20})
+plt.rcParams.update({'font.size': 25})
 '''
 This 
 '''
@@ -25,6 +28,7 @@ class Plot1DHist:
         self.calibration_factor_b1 = 0.652
         self.calibration_factor_b2 = 0.668
         self.silizium = 3.6e-3
+        self.resV = 0.838214/2
         self.b1_mpvs = []
         self.b2_mpvs = []
         self.b1b2_mpvs = []
@@ -32,73 +36,121 @@ class Plot1DHist:
         self.b2_errors = []
         self.b1b2_errors = []
         self.names = []
-    
         self.plot_hist()
         
     def plot_hist(self):
         self.file_names()
+        '''
+        a = amplitude
+        e = position
+        s = width
+        ab = amplitude bg
+        sb = width bg
+        '''
+
+        self.fitting_values = pd.DataFrame({ 'name': [], 'a': [], 'e': [], 's': [], 'a_error': [], 'e_error': [], 's_error': []})
         for file_name in tqdm(self.files, desc="Processing files", unit="file"):
             fig, ax = plt.subplots(figsize=(20, 10))
             data = np.loadtxt(file_name, delimiter=',', comments='#')
             # lese Messzeit aus erste Zeile ein
             with open(file_name, 'r') as file:
                 measurement_time =  int(next(file).split(':')[-1].strip()) / 60 / 60
-            # corrected 
-            B2 = data[:, 2] / measurement_time
-            B1 = data[:, 1] / measurement_time
+                self.measurement_time = measurement_time
+            # B1 and B2 switched because mapping of B1 and B2 is wrong? 
+            B1 = data[:, 2] / measurement_time / self.resV
+            B2 = data[:, 1] / measurement_time / self.resV
             bins = data[:, 0]
-            sum_B1_B2 = data[:, 3] / measurement_time
-            bins_b1 = bins / self.calibration_factor_b1
-            bins_b2 = bins / self.calibration_factor_b2
+            sum_B1_B2 = data[:, 3] / measurement_time / self.resV
+            bins_b1 = bins / self.calibration_factor_b2
+            bins_b2 = bins / self.calibration_factor_b1
             numbers = re.findall(r'\d+', os.path.basename(file_name))
             # den threshold ermitteln
-            if np.sum(B1) > 10:
-                x_b1, y_b1, mpv_b1, e_error_b1 = self.fit_landau(bins_b1,B1)
+            if np.sum(B1) > 3:
+                x_b1, y_b1, mpv_b1, e_error_b1, error, popt = self.fit_landau(bins_b1,B1)# 
+
                 name = fr'$m_{{\text{{TO{numbers[0]}TU{numbers[1]}B1}}}}$ = '
-                ax.step(bins_b1, B1, label= 'B1 mpv: '+name+'('+str(round(mpv_b1,2))+'$\pm$'+str(math.ceil(e_error_b1*100)/100)+')' + ' $keV_{Si}$', where='mid')
-                self.b1_mpvs.append(mpv_b1)
-                self.b1_errors.append(e_error_b1)
+                round_valueb1, round_errorb1 = round_DIN1333(mpv_b1/self.silizium/1000000, e_error_b1/self.silizium/1000000)
+                round_valueb1 = self.multiply_with_error(float(round_valueb1))
+                round_errorb1 = self.multiply_with_error(float(round_errorb1))
+                ax.step(bins_b1 / self.silizium, B1, label= 'B1 mpv: '+name+str(round_valueb1)+'$\pm$'+str(round_errorb1), where='mid')
+                self.b1_mpvs.append(round_valueb1)
+                self.b1_errors.append(round_errorb1)
+                # append all fitting values and errors to dataframe
+                new_row = pd.DataFrame({'name': f'TO{numbers[0]}TU{numbers[1]}B1', 'a': popt[0], 'e': popt[1], 's': popt[2], 'a_error': error[0], 'e_error': error[1], 's_error': error[2]}, index=[0])
+                self.fitting_values = self.fitting_values._append(new_row, ignore_index=True)
 
                 name = fr'$m_{{\text{{TO{numbers[0]}TU{numbers[1]}B2}}}}$ = '
-                x_b2, y_b2, mpv_b2, e_error_b2 = self.fit_landau(bins_b2,B2)
-                ax.step(bins_b2, B2, label= 'B2 mpv: '+name+'('+str(round(mpv_b2,2))+'$\pm$'+str(math.ceil(e_error_b2*100)/100)+')' + ' $keV_{Si}$' , where='mid')
-                self.b2_mpvs.append(mpv_b2)
-                self.b2_errors.append(e_error_b2)
+                x_b2, y_b2, mpv_b2, e_error_b2, error, popt = self.fit_landau(bins_b2,B2)
+                round_valueb2, round_errorb2 = round_DIN1333(mpv_b2/self.silizium/1000000, e_error_b2/self.silizium/1000000)
+                round_valueb2 = self.multiply_with_error(float(round_valueb2))
+                round_errorb2 = self.multiply_with_error(float(round_errorb2))
+                ax.step(bins_b2/ self.silizium, B2, label= 'B2 mpv: '+name+str(round_valueb2)+'$\pm$'+str(round_errorb2) , where='mid')
+                self.b2_mpvs.append(round_valueb2)
+                self.b2_errors.append(round_errorb2)
+                new_row = pd.DataFrame({'name': f'TO{numbers[0]}TU{numbers[1]}B2', 'a': popt[0], 'e': popt[1], 's': popt[2], 'a_error': error[0], 'e_error': error[1], 's_error': error[2]}, index=[0])
+                self.fitting_values = self.fitting_values._append(new_row, ignore_index=True)
 
                 name = fr'$m_{{\text{{TO{numbers[0]}TU{numbers[1]}B1B2}}}}$ = '
-                x_b1b2, y_b1b2, mpv_b1b2, e_error_b1b2 = self.fit_landau(bins,sum_B1_B2, sum = True)
-                ax.step(bins, sum_B1_B2, label= 'B1+B2 mpv: '+name+'('+str(round(mpv_b1b2,2))+'$\pm$'+str(math.ceil(e_error_b1b2*100)/100)+')' + ' $keV_{Si}$' , where='mid')
-                self.b1b2_mpvs.append(mpv_b1b2)
-                self.b1b2_errors.append(e_error_b1b2)
+                x_b1b2, y_b1b2, mpv_b1b2, e_error_b1b2, error, popt = self.fit_landau(bins,sum_B1_B2, sum = True)
+                round_valueb1b2, round_errorb1b2 = round_DIN1333(mpv_b1b2/self.silizium/1000000, e_error_b1b2/self.silizium/1000000)
+                round_valueb1b2 = self.multiply_with_error(float(round_valueb1b2))
+                round_errorb1b2 = self.multiply_with_error(float(round_errorb1b2))
+                ax.step(bins/ self.silizium, sum_B1_B2, label= 'B1+B2 mpv: '+name+str(round_valueb1b2)+'$\pm$'+str(round_errorb1b2)  , where='mid')
+                self.b1b2_mpvs.append(round_valueb1b2)
+                self.b1b2_errors.append(round_errorb1b2)
                 self.names.append(os.path.basename(file_name))
+                new_row = pd.DataFrame({'name': f'TO{numbers[0]}TU{numbers[1]}B1B2', 'a': popt[0], 'e': popt[1], 's': popt[2], 'a_error': error[0], 'e_error': error[1], 's_error': error[2],}, index = [0])
+                self.fitting_values = self.fitting_values._append(new_row, ignore_index=True)
 
-                ax.plot(x_b1, y_b1, linestyle='--', color='black')
-                ax.plot(x_b2, y_b2, linestyle='--', color='black')
-                ax.plot(x_b1b2, y_b1b2, linestyle='--', color='black')
+                ax.plot(x_b1/ self.silizium, y_b1, linestyle='--', color='black')
+                ax.plot(x_b2/ self.silizium, y_b2, linestyle='--', color='black')
+                ax.plot(x_b1b2/ self.silizium, y_b1b2, linestyle='--', color='black')
             else:
-                ax.step(bins_b1, B1, label='B1', where='mid')
-                ax.step(bins_b2, B2, label='B2', where='mid',)
-                ax.step(bins, sum_B1_B2, label='B1+B2', where='mid')
+                ax.step(bins_b1/ self.silizium, B1, label='B1', where='mid')
+                ax.step(bins_b2/ self.silizium, B2, label='B2', where='mid',)
+                ax.step(bins/ self.silizium, sum_B1_B2, label='B1+B2', where='mid')
 
             ax.set_yscale('log')
-            ax.set_xlim(0, 500)
-            ax.set_ylim(bottom = 0.9/ measurement_time)
-            plt.title('Histogramm muon coincidence of B1, B2, TO' + str(numbers[0]) + ' and ' + 'TU' + str(numbers[1]))
+            ax.set_xlim(0, 500/ self.silizium)
+            ax.set_ylim(bottom = 0.9/ measurement_time/ self.resV)
+            # set upper limit of y
+            ax.set_ylim(top = 1)
+            plt.title('1D histogram coincidence of TO' + str(numbers[0]) + ' and ' + 'TU' + str(numbers[1]))
             plt.legend()
-            plt.xlabel(r'Energy [$keV_{Si}$]')
+            plt.xlabel('emitted photoelectrons')
             plt.grid()
-            plt.ylabel('Counts [1/h/mV]')
+            plt.ylabel('counts [1/h/mV]')
             # Load and display the image in the top-left corner
             image_path = 'DiodenBezeichnung.png'
             setup = image.imread(image_path)
-            imagebox = OffsetImage(setup, zoom=0.3)
-            ab = AnnotationBbox(imagebox, (0.8, 0.55), xycoords='figure fraction', frameon=False)
+            imagebox = OffsetImage(setup, zoom=0.4)
+            ab = AnnotationBbox(imagebox, (0.8, 0.45), xycoords='figure fraction', frameon=False)
             fig.add_artist(ab)
 
             file_path = self.target_directory_muon+'plot'+'/'+os.path.basename(file_name[:-4]) + '.pdf'
             plt.savefig(file_path)
+            #plt.show()
             plt.close()
         self.save_txt_mpvs()
+        self.save_fitting_values()
+    # Function to round float values to 3 digits
+    def round_float_values(self, value):
+        if isinstance(value, float):
+            return round(value, 3)
+        return value
+    def multiply_with_error(self,value):
+        # Set precision to a sufficiently high value
+        getcontext().prec = 50  # You can adjust the precision as needed
+        # Convert values to Decimal
+        value_decimal = Decimal(str(value))
+        factor_decimal = Decimal(str(1000000))
+        # Perform multiplication with arbitrary precision
+        new_value = value_decimal * factor_decimal
+        return float(str(new_value).rstrip('0').rstrip('.'))
+    def save_fitting_values(self):
+        # rround like in the upper code all values corresponding to their errors
+
+        self.fitting_values.to_csv(self.target_directory_muon+'/..'+'/'+'fitting_values.txt', index=False)
 
     def save_txt_mpvs(self):
         # create pandas dataframe
@@ -112,8 +164,9 @@ class Plot1DHist:
         self.files = file_names_muon
 
     def fit_landau(self, x_data, y_data, sum = False):
-        mid_cut_low = 20
-        mid_cut_high = 200
+        mid_cut_low = 30
+        # old cut 130
+        mid_cut_high = 130
         if sum:
             mid_cut_low = 80
             mid_cut_high = 400
@@ -122,39 +175,48 @@ class Plot1DHist:
         y_s = []
         mpv_s = []
         e_errors = []
-        for lower_cut in range(0,10,2):
-            for higher_cut in range(0,50,10):
+        error_list = []
+        popt_list = []
+        for lower_cut in range(0,20,2):
+            for higher_cut in range(0,50,2):
                 try:   
                     mask = (x_data > mid_cut_low+lower_cut) & (x_data < mid_cut_high+higher_cut)
-                    new_x, y, mpv, e_error = self.fit_function_landau(x_data[mask], y_data[mask])
-                    error = np.sum(np.abs(y_data[mask][1:] - y))
+                    #print(x_data[mask])
+                    new_x, y, mpv, e_error, errors_all, popt_all = self.fit_function_landau(x_data[mask], y_data[mask])
+                    #error = np.sum(np.abs(y_data[mask][1:] - y))
+                    error_list.append(errors_all)
+                    popt_list.append(popt_all)
                     x_s.append(new_x)
                     y_s.append(y)
                     mpv_s.append(mpv)
-                    errors.append(error)
+                    errors.append(e_error)
                     e_errors.append(e_error)
                 except Exception:
                     pass
         # find best fit
         best_fit = np.argmin(errors)
-        return x_s[best_fit], y_s[best_fit], mpv_s[best_fit], e_errors[best_fit]
+        #print(e_errors)
+        return x_s[best_fit], y_s[best_fit], mpv_s[best_fit], e_errors[best_fit], error_list[best_fit], popt_list[best_fit]
 
     def fit_function_landau(self, x_data, y_data):
         # fit on middle of the bins and not on the edges
-        x_data = x_data[:-1] + (x_data[1] - x_data[0]) / 2
+        x_data = x_data[:-1]
         y_data = y_data[1:]        
         # First iteration
-        sigma = np.sqrt(y_data+1)
-        p0 = [np.max(y_data), x_data[np.argmax(y_data)]-10, 10, np.max(y_data), 100]
-        popt, pcov = curve_fit(self.lan, x_data, y_data, p0=p0, maxfev = 1000, sigma=sigma, absolute_sigma=True)
-
-        a_error, e_error, s_error, ab_error, sb_error = np.sqrt(np.diag(pcov))
+        sigma = np.sqrt(y_data+1/self.measurement_time/self.resV)
+        p0 = [np.max(y_data), x_data[np.argmax(y_data)], 10]
+        #print(np.max(y_data))
+        popt, pcov = curve_fit(self.lan, x_data, y_data, p0=p0, maxfev = 2000, sigma=sigma, absolute_sigma=True)
+        a_error, e_error, s_error = np.sqrt(np.diag(pcov))
         #new_x = np.linspace(x_data[0], x_data[-1], 100)
         y = self.lan(x_data, *popt)
+        # all errors from covariance matrix
+        error = np.sqrt(np.diag(pcov))
+
         mpv = popt[1]
         #mpv = np.argmax(y)
         #mpv = new_x[mpv]
-        return x_data, y, mpv, e_error
+        return x_data, y, mpv, e_error, error, popt
     
     def mips(self,x, a,e,s):
         return a*self.landau((x-e)/s)
@@ -162,11 +224,7 @@ class Plot1DHist:
     def landau(self,l):
         return np.sqrt(np.exp(-l-np.exp(-l))/2/np.pi)
 
-    def bg(self,x, a,e,s):
-        return a*np.exp((e-x)/s)
-
-    
-    def lan(self,x, a,e,s,ab,sb):
+    def lan(self,x, a,e,s):
         '''
         a = amplitude
         e = position
@@ -174,7 +232,8 @@ class Plot1DHist:
         ab = amplitude bg
         sb = width bg
         '''
-        return self.mips(x, a,e,s) + self.bg(x, ab, e, sb)
+        return self.mips(x, a,e,s) 
+    
     
 # Create an instance of the class
 source_file = '2023-12-21-AHBGOC-14-langzeit.dat'

plot_1d_histogramm_xray.py

@@ -6,14 +6,24 @@ from scipy.optimize import curve_fit
 from matplotlib.ticker import MultipleLocator
 #ignore warnings
 import warnings
+import math
+from decimal import Decimal, getcontext
+import pandas as pd
+from round_DIN1333 import round_DIN1333
 from toms_fit_function import compton
 from scipy.constants import m_e,c, elementary_charge
 warnings.filterwarnings("ignore")
 plt.rcParams.update({'font.size': 20})
 
 class Plot1DHist:
-    def __init__(self, dat_name = '2023-12-06-AHBGOC-Bi207-07.dat'):
+    def __init__(self, dat_name = '2023-12-01-AHBGOC-Bi207-06.dat'):
+        # 6er ist oben
+        # 7er ist unten
+        #dat_name = '2023-12-06-AHBGOC-Bi207-07.dat'
+        #dat_name = '2023-12-11-AHBGOC-Bi207-08.dat'
         self.target_directory_muon = f'/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results{dat_name[:-4]}/1d_hist_xray'
+        # unten 
+
         os.makedirs(self.target_directory_muon+'plot', exist_ok=True)
         # energy first compton edge in eV
         self.energy_first_compton_edge = 569698
@@ -23,6 +33,11 @@ class Plot1DHist:
 
     def plot_hist(self):
         self.file_names()
+        fitting_values = []
+        fitting_errors = []
+        fitting_names = []
+        self.fitting_values = pd.DataFrame({'name': [], 'a': [],'b':[], 'e':[], 's':[], 'r':[], 'a_error': [], 'e_error':[], 's_error':[], 'r_error':[]})
+        fitting_values = []
         fig, ax = plt.subplots(figsize=(20, 10))
         for file_name in tqdm(self.files, desc="Processing files", unit="file"):
             data = np.loadtxt(file_name, delimiter=',', skiprows=2)
@@ -31,42 +46,44 @@ class Plot1DHist:
                 measurement_time =  int(next(file).split(':')[-1].strip()) / 60 / 60
             B1 = data[:, 1] / measurement_time
             bins = data[:, 0]
-            ax.step(bins, B1, where='mid', label=os.path.basename(file_name[:-9]))
+            #ax.step(bins, B1, where='mid', label=os.path.basename(file_name[:-9]))
 
-            # no gaussian fit because it is a big pain in the ass
-            '''# Gaussian fit
-            if 'B' in os.path.basename(file_name):
-                new_x, y,u, u_error = self.fit_gauss(bins,B1, BGO=True)
-                ax.plot(new_x, y, linestyle='--', color='black', label='Gaussian fit: '+str(round(u,2))+'$\pm$'+str(round(u_error,2))+'mV/keV')
-            else:
-                new_x, y,u, u_error, u1, u1_error = self.fit_gauss(bins,B1)
-                ax.plot(new_x, y, linestyle='--', color='black', label='Gaussian fit: '+str(round(u,2))+'$\pm$'+str(round(u_error,2))+' mV/keV'+'\n'+str(round(u1,2))+'$\pm$'+str(round(u1_error,2))+' mV/keV')
-            '''
             # Compton edge
-            x,y, e, e_error = self.fit_compton(bins,B1)
-            ax.plot(x,y, label=rf'Fit Compton edge: $u_{{\text{{e,bgo{os.path.basename(file_name)[1]}}}}}$'+' = ('+str(round(self.calculate_calibration_first_edge(e),3))+'$\pm$'+str(round(self.calculate_calibration_first_edge(e_error),3))+') mV/keV', color='black', linestyle='--')
-            print('basename: ', os.path.basename(file_name))
-
+            x,y, e, e_error,best_vals, errors = self.fit_compton(bins,B1)
+            fitting_values.append(best_vals)
+            fitting_errors.append(errors)
+            fitting_names.append(file_name)
+            print('edge calibration: ', self.calculate_calibration_first_edge(e), 'e_error: ', self.calculate_calibration_first_edge(e_error))
+            rounded_edge, rounded_edge_error = round_DIN1333(self.calculate_calibration_first_edge(e), self.calculate_calibration_first_edge(e_error))
+            ax.step(bins, B1, where='mid', label=os.path.basename(file_name[:-9]) + rf': fit compton edge: $u_{{\text{{edge}}}}$'+' = ('+str(rounded_edge)+'$\pm$'+str(rounded_edge_error)+') mV/keV')
+            ax.plot(x,y, color='black', linestyle='--')
+            a,b,e,s,r = np.round(best_vals, 3)
+            a_error ,b_error, e_error, s_error, r_error = np.round(errors, 3)
+            new_row = pd.DataFrame({'name': os.path.basename(file_name[:-9]), 'a': a,'b': b, 'e': round(e,4), 's': s, 'r': r, 'a_error': a_error, 'b_error': b_error,'e_error': round(e_error,4), 's_error': s_error, 'r_error': r_error}, index = [0])
+            #new_row = pd.DataFrame({'name': os.path.basename(file_name[:-9]), 'e': round(self.calculate_calibration_first_edge(e),3)}, index = [0])
+            self.fitting_values = self.fitting_values._append(new_row, ignore_index = True)
         ax.set_yscale('log')
         ax.set_xlim(0, 800)
         ax.set_ylim(bottom = 0.9/ measurement_time)
         plt.title('Calibration from run '+ self.dat_name)
         plt.legend()
-        plt.xlabel('A[mV]')
+        plt.xlabel('A [mV]')
         ax.xaxis.set_minor_locator(MultipleLocator(10))
         plt.grid()
         plt.ylabel('counts [1/h/mV]')
         file_path = self.target_directory_muon+'plot'+'/'+os.path.basename(file_name[:-4]) + '.pdf'
         print('file_path: ', file_path)
-        plt.savefig(file_path)
+        plt.savefig('fitting_values_6_B1.pdf')
         plt.show()
         plt.close()
+        self.fitting_values.to_csv('fitting_values_unten_all.csv')
 
     def file_names(self):
         file_names_muon = os.listdir(self.target_directory_muon)
         # remove files which contain an 'O'
-        file_names_muon = [file_name for file_name in file_names_muon if 'B' in file_name]
+        file_names_muon = [file_name for file_name in file_names_muon if 'U232' in file_name or 'B' in file_name]
         file_names_muon = [self.target_directory_muon + '/' + file_name for file_name in file_names_muon]
+        print('file_names_muon: ', file_names_muon)
         self.files = file_names_muon
     
     def calculate_calibration_first_edge(self,e):
@@ -79,8 +96,8 @@ class Plot1DHist:
         pass
     
     def fit_compton(self,x,y):
-        mid_cut_low = 200
-        mid_cut_high = 310
+        mid_cut_low = 220
+        mid_cut_high = 300
         init_vals = [3, 1, 260, 50, 100]
 
         mask = (x > mid_cut_low) & (x < mid_cut_high)
@@ -90,7 +107,17 @@ class Plot1DHist:
         error_e  = errors[1]
         x = np.linspace(x[0]-100,x[-1]+100,10000)
         y = compton(x, *best_vals)
-        return x,y, best_vals[2], error_e
+        return x,y, best_vals[2], error_e,best_vals, errors   
+    
+    def multiply_with_error(self,value):
+        # Set precision to a sufficiently high value
+        getcontext().prec = 50  # You can adjust the precision as needed
+        # Convert values to Decimal
+        value_decimal = Decimal(str(value))
+        factor_decimal = Decimal(str(100000))
+        # Perform multiplication with arbitrary precision
+        new_value = value_decimal * factor_decimal
+        return float(str(new_value).rstrip('0').rstrip('.'))
         
     def fit_gauss(self,x,y, init_vals = [10, 0.66, 1, 4,8,10,10,10,50,20], BGO = False):
         #bezeichnung                     s,u,a_alpha,a_beta, m,e ,r, s1, u1, a1
@@ -105,8 +132,6 @@ class Plot1DHist:
                     if BGO:
                         best_vals, covar = curve_fit(self.hole_gauss_BGO, xdata=x[mask], ydata=y[mask], p0=init_vals[:-3][:3] + init_vals[:-3][4:], maxfev = 5000)
                         errors.append(np.sum((y[mask] - self.hole_gauss_BGO(x[mask], *best_vals))**2))
-
-                        print('best_vals', best_vals)
                     else:
                         best_vals, covar = curve_fit(self.hole_gauss, xdata=x[mask], ydata=y[mask], p0=init_vals, maxfev = 5000)
                         errors.append(np.sum((y[mask] - self.hole_gauss(x[mask], *best_vals))**2))
@@ -115,10 +140,7 @@ class Plot1DHist:
                     uppers.append(upper_limit)  
                 except Exception as e:
                     pass
-                    '''print(f"Error type: {type(e).__name__}")
-                    print(f"Error message: {e}")
-                    import traceback
-                    traceback.print_exc()'''
+
 
         # find best fit
         best_fit = np.argmin(np.array(errors))  
@@ -180,4 +202,4 @@ class Plot1DHist:
     
 # Create an instance of the class
 source_file = '2023-12-11-AHBGOC-Bi207-08.dat'
-plot_instance = Plot1DHist(dat_name = source_file)
+plot_instance = Plot1DHist()

plot_2d_histogramm_muon.py

@@ -4,11 +4,13 @@ import numpy as np
 from tqdm import tqdm
 from matplotlib import cm
 import re
+import pandas as pd
 from scipy.optimize import curve_fit
 import matplotlib as mpl
+from decimal import Decimal, getcontext
 import matplotlib.image as image
 from matplotlib.offsetbox import (OffsetImage, AnnotationBbox)
-
+from round_DIN1333 import  round_DIN1333
 # Set a larger default font size
 mpl.rcParams['font.size'] = 25
 
@@ -19,8 +21,10 @@ class Plot2DHist:
                 self.path_to_file = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/2dhistmuons/'
                 os.makedirs(self.target_directory+'plot', exist_ok=True)
                 self.minV = -100
-                self.calibration_factor_b1 = 0.652
-                self.calibration_factor_b2 = 0.668
+                # b1 and b2 seem to be correct
+                self.calibration_factor_b1 = 0.6524
+                self.calibration_factor_b2 = 0.6684
+                self.silizium = 3.6e-3
                 self.maxV = 3500
                 self.resV = 0.838214/2
                 self.bins = np.arange(self.minV, self.maxV, self.resV)
@@ -31,59 +35,84 @@ class Plot2DHist:
         def create_file_names(self):
                 self.file_names = os.listdir(self.target_directory)
 
+        def multiply_with_error(self,value, multiplier):
+                # Set precision to a sufficiently high value
+                getcontext().prec = 50  # You can adjust the precision as needed
+                # Convert values to Decimal
+                value_decimal = Decimal(str(value))
+                factor_decimal = Decimal(str(multiplier))
+                # Perform multiplication with arbitrary precision
+                new_value = value_decimal * factor_decimal
+                return float(str(new_value).rstrip('0').rstrip('.'))
+        
         def plot_hist(self):
+                fiiting_values = pd.DataFrame(columns=['Coincidence', 'slope', 'intercept in measured charge Q', 'slope_error', 'intercept_error in measured charge Q'])
                 for i in tqdm(self.file_names):
                         numbers = re.findall(r'\d+', i)
-                        fig, ax = plt.subplots(figsize=(20, 10))
-                        self.data = np.loadtxt(self.path_to_file + i, delimiter=',')
-                        xlim = (0,250)
-                        ylim = (0,250)
-                        cmap = cm.get_cmap('jet', 256)
-                        cmap.set_under('white')
+                        bedingung = numbers[0] == '9' and numbers[1] == '9'
+                        if True:
+                                fig, ax = plt.subplots(figsize=(20, 10))
+                                self.data = np.loadtxt(self.path_to_file + i, delimiter=',')
+                                xlim = (0,250)
+                                ylim = (0,250)
+                                cmap = cm.get_cmap('jet', 256)
+                                cmap.set_under('white')
                                 
-                        self.bins_new = self.bins[(self.bins >= xlim[0] ) & (self.bins <= xlim[1]) & (self.bins >= ylim[0]) & (self.bins <= ylim[1])]
-                        self.bins_new_b1 = self.bins_new / self.calibration_factor_b1
-                        self.bins_new_b2 = self.bins_new / self.calibration_factor_b2
+                                self.bins_new = self.bins[(self.bins >= xlim[0] ) & (self.bins <= xlim[1]) & (self.bins >= ylim[0]) & (self.bins <= ylim[1])]
+                                self.bins_new_b1 = self.bins_new / self.calibration_factor_b2
+                                self.bins_new_b2 = self.bins_new / self.calibration_factor_b1
                                 
-                        self.data = self.data[(self.bins[:-1] >= xlim[0]) & (self.bins[:-1] <= xlim[1]), :]
-                        self.data = self.data[:, (self.bins[:-1] >= ylim[0]) & (self.bins[:-1] <= ylim[1])].T
+                                self.data = self.data[(self.bins[:-1] >= xlim[0]) & (self.bins[:-1] <= xlim[1]), :]
+                                # y axis corresponds to b1 and x-axis corresponds to b2
+                                self.data = self.data[:, (self.bins[:-1] >= ylim[0]) & (self.bins[:-1] <= ylim[1])]
 
-                        # ist vertauscht weil die matrix so gespeichert wird. Würde man self.data Transponieren, wäre es richtig               
-                        im = ax.pcolormesh(self.bins_new_b1, self.bins_new_b2, self.data, cmap=cmap, shading='auto', vmin=0.6)
-                        x = y = np.arange(xlim[0],xlim[1], 1)
-                        ax.plot(x, y, color='black',linewidth=2, linestyle='--',label='line with slope 1 and intercept 0')
-                        # fit line
-                        if np.sum(self.data) > 500:
-                                self.line_data_from_matrice()
-                                popt, pcov = curve_fit(self.line, self.x, self.y,p0=[1,0], method='lm')
-                                a_error = np.sqrt(pcov[0][0])
-                                b_error = np.sqrt(pcov[1][1])
-                                ax.plot(x, self.line(x, *popt), color='red', linestyle='--', linewidth=2, label=f'line fit:\nslope s = ({round(popt[0],3)} ± {round(a_error,3)})\nintercept b = ({round(popt[1],1)} ± {round(b_error,1)})'+' $keV_{Si}$')
+                                # ist vertauscht weil die matrix so gespeichert wird. Würde man self.data Transponieren, wäre es richtig               
+                                im = ax.pcolormesh(self.bins_new_b2/ self.silizium, self.bins_new_b1/ self.silizium, self.data, cmap=cmap, shading='auto', vmin=0.6)
+                                x = y = np.arange(xlim[0],xlim[1], 1)
+                                ax.plot(x/ self.silizium, y/ self.silizium, color='black',linewidth=2, linestyle='--',label='line with slope 1 and intercept 0')
+                                # fit line
+                                if np.sum(self.data) > 500:
+                                        name = 'TO' + numbers[0] + '_TU' + numbers[1]
+                                        # create pandas dataframe with a,b, a and b errors and save as csv
+                                        self.line_data_from_matrice()
+                                        popt, pcov = curve_fit(self.line, self.x, self.y,p0=[1,0])
+                                        a_error = np.sqrt(pcov[0][0])
+                                        b_error = np.sqrt(pcov[1][1])
+                                        round_a_value, round_a_error = round_DIN1333(popt[0]/10, a_error/10) 
+                                        round_b_value, round_b_error = round_DIN1333(popt[1]/self.silizium/10000, b_error/self.silizium/10000)
+                                        round_a_value = self.multiply_with_error(round_a_value,10)
+                                        round_b_value = self.multiply_with_error(round_b_value,10000)
+                                        round_a_error = self.multiply_with_error(round_a_error,10)
+                                        round_b_error = self.multiply_with_error(round_b_error,10000)   
+                                        ax.plot(x / self.silizium, self.line(x, *popt)/ self.silizium, color='red', linestyle='--', linewidth=2, label=f'line fit:\nslope s = {round_a_value} ± {round_a_error}\nintercept b = {round_b_value} ± {round_b_error}')
+                                        new_row = {'Coincidence': name, 'slope': round_a_value, 'intercept in measured charge Q': round_b_value, 'slope_error': round_a_error, 'intercept_error in measured charge Q': round_b_error}
+                                        fiiting_values = fiiting_values._append(new_row, ignore_index=True)
+                                ax.set_xlim(xlim[0]/ self.silizium, xlim[1]/ self.silizium)
+                                ax.set_ylim(ylim[0]/ self.silizium, ylim[1]/ self.silizium)
+                                plt.grid()
+                                plt.title(os.path.basename(self.path_to_file))
+                                plt.xlabel('B2 emitted photoelectrons')
+                                plt.ylabel('B1 emitted photoelectrons')
+                                plt.legend(loc = 'upper left')
+                                plt.title(f'2D histogram coincidence of TO{numbers[0]} and TU{numbers[1]}')
+                                # Create colorbar using the image object
+                                cbar = plt.colorbar(im)
+                                cbar.set_label('Counts')  # Set your colorbar label
+                                coinc = re.findall(r'\d+', i)
+                                file_name = f'TO{coinc[0]}_TU{coinc[1]}_2dhist.png'
 
-                        ax.set_xlim(xlim[0], xlim[1])
-                        ax.set_ylim(ylim[0], ylim[1])
-                        plt.grid()
-                        plt.title(os.path.basename(self.path_to_file))
-                        plt.xlabel('B1 energy deposit in [$keV_{Si}$]')
-                        plt.ylabel('B2 energy deposit in [$keV_{Si}$]')
-                        plt.legend(loc = 'upper left')
-                        plt.title(f'2D histogramm coincidence of TO{numbers[0]}, TU{numbers[1]}, B1 and B2')
-                        # Create colorbar using the image object
-                        cbar = plt.colorbar(im)
-                        cbar.set_label('Counts')  # Set your colorbar label
-                        coinc = re.findall(r'\d+', i)
-                        file_name = f'TO{coinc[0]}_TU{coinc[1]}_2dhist.png'
-
-                        # add setup 
-                        image_path = 'DiodenBezeichnung.png'
-                        setup = image.imread(image_path)
-                        imagebox = OffsetImage(setup, zoom=0.3)
-                        ab = AnnotationBbox(imagebox, (0.65, 0.3), xycoords='figure fraction', frameon=False)
-                        fig.add_artist(ab)
-
-                        plt.savefig(self.target_directory+'plot'+'/'+file_name)
-                        # clear memory
-                        plt.close('all')
+                                # add setup 
+                                image_path = 'DiodenBezeichnung.png'
+                                setup = image.imread(image_path)
+                                imagebox = OffsetImage(setup, zoom=0.3)
+                                ab = AnnotationBbox(imagebox, (0.65, 0.3), xycoords='figure fraction', frameon=False)
+                                fig.add_artist(ab)
+                                #print('path', self.target_directory+'plot'+'/'+file_name)
+                                plt.savefig(self.target_directory+'plot'+'/'+file_name)
+                                #plt.show()
+                                # clear memory
+                                plt.close('all')
+                fiiting_values.to_csv('fitting_values_2dHistogramm.csv', index=False)
 
         def line(self,x,a,b):
                 return a*x+b
@@ -93,16 +122,13 @@ class Plot2DHist:
                 # iterate over each entry in the matrix and append the x and y value to a list, so the corresponding bins will add
                 x = []
                 y = []
-                # Calculate the total number of iterations
                 total_iterations = int(np.sum(self.data))
                 print('total counts: ', total_iterations)
-                # Use tqdm for a progress bar
                 for index_row, row in enumerate(self.data):
                         for index_column, column in enumerate(row):
                                 for times in range(int(self.data[index_row][index_column])):
-                                        # x is B1 so 
-                                        x.append(self.bins_new_b1[index_column])
-                                        y.append(self.bins_new_b2[index_row])
+                                        x.append(self.bins_new_b2[index_column])
+                                        y.append(self.bins_new_b1[index_row])
 
                 self.x = np.array(x)
                 self.y = np.array(y)

plot_after_calibration.py

@@ -0,0 +1,116 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import os
+from scipy.optimize import curve_fit
+import math
+plt.rcParams.update({'font.size': 25})
+
+class AfterCalibration:
+    def __init__(self) -> None:
+        self.target_path = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/1d_hist_muon_no_coinc'
+        self.calibrationb1 = 0.652
+        self.calibrationb2 = 0.668
+        self.silizium = 3.6e-3
+        self.fet_file_names()
+        self.plot_hist()
+
+    def fet_file_names(self):
+        self.file_names = os.listdir(self.target_path)
+
+    def plot_hist(self):
+        plt.figure(figsize=(20,10))
+        for file_name in self.file_names:
+            if 'B1' in os.path.basename(file_name):
+                # fit landau
+                data = np.loadtxt(self.target_path + '/' + file_name, delimiter=',', skiprows=2, comments='#')
+                bins = data[:, 0] / self.calibrationb1
+                hist = data[:, 1]
+                x,y,mpv,e_error = self.fit_landau(bins, hist)
+                plt.step(bins/self.silizium, hist, label='B1 mpv: $m_{B1}$ = '+f'{round(mpv/self.silizium)} $\pm$ {math.ceil(e_error/self.silizium)}')
+                plt.plot(x/self.silizium,y, linestyle='--', color='black')
+            if 'B2' in os.path.basename(file_name):
+                data = np.loadtxt(self.target_path + '/' + file_name, delimiter=',', skiprows=2, comments='#')
+                bins = data[:, 0] / self.calibrationb2
+                hist = data[:, 1]
+                x,y,mpv,e_error = self.fit_landau(bins, hist)
+                plt.step(bins/self.silizium, hist, label='B2 mpv: $m_{B2}$ = '+f'{round(mpv/self.silizium)} $\pm$ {math.ceil(e_error/self.silizium)} ')
+                plt.plot(x/self.silizium,y, linestyle='--', color='black')
+
+            plt.xlabel('emitted photoelectrons')
+            plt.ylabel('Counts')
+            plt.xlim(30/self.silizium, 150/self.silizium)
+            plt.ylim(5000,30000)
+            plt.legend()
+            plt.title('Histogram muon 2023-12-21-AHBGOC-14-langzeit.dat')
+            plt.yscale('log')
+        plt.grid()
+        plt.savefig(self.target_path + '/../' + file_name[:-4] + '.pdf')
+        plt.show()
+    def fit_landau(self, x_data, y_data, sum = False):
+        mid_cut_low = 40
+        mid_cut_high = 250
+        if sum:
+            mid_cut_low = 80
+            mid_cut_high = 400
+        errors = []
+        x_s = []
+        y_s = []
+        mpv_s = []
+        e_errors = []
+        for lower_cut in range(0,10,2):
+            for higher_cut in range(0,100,10):
+                try:   
+                    mask = (x_data > mid_cut_low+lower_cut) & (x_data < mid_cut_high+higher_cut)
+                    new_x, y, mpv, e_error = self.fit_function_landau(x_data[mask], y_data[mask])
+                    error = np.sum(np.abs(y_data[mask][1:] - y))
+                    x_s.append(new_x)
+                    y_s.append(y)
+                    mpv_s.append(mpv)
+                    errors.append(error)
+                    e_errors.append(e_error)
+                except Exception:
+                    pass
+        # find best fit
+        best_fit = np.argmin(errors)
+        return x_s[best_fit], y_s[best_fit], mpv_s[best_fit], e_errors[best_fit]
+
+    def fit_function_landau(self, x_data, y_data):
+        # fit on middle of the bins and not on the edges
+        x_data = x_data[:-1] + (x_data[1] - x_data[0]) / 2
+        y_data = y_data[1:]        
+        # First iteration
+        sigma = np.sqrt(y_data+1)
+        p0 = [np.max(y_data), x_data[np.argmax(y_data)]-10, 10, np.max(y_data), 100]
+        popt, pcov = curve_fit(self.lan, x_data, y_data, p0=p0, maxfev = 1000, sigma=sigma, absolute_sigma=True)
+        
+        a_error, e_error, s_error, ab_error, sb_error = np.sqrt(np.diag(pcov))
+        #new_x = np.linspace(x_data[0], x_data[-1], 100)
+        y = self.lan(x_data, *popt)
+        mpv = popt[1]
+        #mpv = np.argmax(y)
+        #mpv = new_x[mpv]
+        return x_data, y, mpv, e_error
+    
+    def mips(self,x, a,e,s):
+        return a*self.landau((x-e)/s)
+    
+    def landau(self,l):
+        return np.sqrt(np.exp(-l-np.exp(-l))/2/np.pi)
+
+    def bg(self,x, a,e,s):
+        return a*np.exp((e-x)/s)
+
+    
+    def lan(self,x, a,e,s,ab,sb):
+        '''
+        a = amplitude
+        e = position
+        s = width
+        ab = amplitude bg
+        sb = width bg
+        '''
+        return self.mips(x, a,e,s) + self.bg(x, ab, e, sb)
+
+
+AfterCalibration()

plot_check_random_hits.py

@@ -0,0 +1,38 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import os
+plt.rcParams.update({'font.size': 18})
+class CheckYield:
+    def __init__(self):
+        self.file_path_only_trigger = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/1d_hist_muon_only_trigger_diodes/'
+        self.file_path_all_trigger = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/1d_hist_muon/'
+        self.file_names = os.listdir(self.file_path_only_trigger)
+        self.chech_ratio_counts_and_plot()
+
+    def chech_ratio_counts_and_plot(self):
+        ratios = []
+        labels = []
+        plt.figure(figsize=(20, 10))
+        for i in self.file_names:
+            labels.append(i[:-5])
+            data_only_trigger = np.loadtxt(self.file_path_only_trigger + i, comments='#', delimiter=',')
+            data_all_trigger = np.loadtxt(self.file_path_all_trigger + i, comments='#', delimiter=',')
+            hist_only_trigger = data_only_trigger[:,1]
+            hist_all_trigger = data_all_trigger[:,1]
+            ratios.append(abs((np.sum(hist_only_trigger)/np.sum(hist_all_trigger)-1)*100))
+            # bar plot of ratios
+        # how much percent higher
+        print('Overall percent higher: ', np.mean(ratios))
+        # sort labels and ratios alphabetically
+        labels, ratios = zip(*sorted(zip(labels, ratios)))
+        plt.bar(labels, ratios)
+        plt.xticks(rotation=90)
+        plt.title('Ratio of counts in 1d histogram with only trigger diodes and trigger diodes plus BGO diodes')
+        plt.ylabel('Ratio of counts')
+        plt.yscale('log')
+        plt.grid()
+        plt.savefig('ratio_counts.pdf', bbox_inches='tight')
+        plt.show()
+
+CheckYield()

plot_hist.py

@@ -8,6 +8,8 @@ from scipy.optimize import curve_fit
 from toms_fit_function import compton
 import plotly.express as px
 import plotly.graph_objects as go
+from itertools import cycle
+plt.rcParams.update({'font.size': 25})
 
 '''
 data which i created with irena:
@@ -25,19 +27,55 @@ data which i created with irena:
 -rw-r--r-- 1 athena  athena     20742656 Dec 21 16:14 2023-12-21-AHBGOC-12.dat -> testmessung
 -rw-r--r-- 1 athena  athena     13955200 Dec 21 17:36 2023-12-21-AHBGOC-13.dat -> testmessung
 -rw-r--r-- 1 athena  athena    122962368 Dec 26 08:40 2023-12-21-AHBGOC-14-langzeit.dat -> Bachelorarbeitsmessung BGO 2 wurde benutzt
+2023-01-16-AHBGOC-rau-15.dat
 '''
 '''
 This code downloads the .dat from asterix and plots the histogramm. to run the code on your own machine, some changes need to be done.
 '''
 
 class Hist:
-    def __init__(self, delete = True, dat_name = '2023-12-21-AHBGOC-14-langzeit.dat'):
+    def __init__(self, delete = False, dat_name = '2023-12-21-AHBGOC-14-langzeit.dat'):
         self.rem_path = '/data/asterix/athena/arm/irena/ahbgo/'
+        #dat_name = '2023-01-16-AHBGOC-rau-15.dat'
+        #dat_name = '2023-01-19-AHBGOC-rau-16.dat'
         self.dat_name = dat_name
         # kalibrationfactor in mV/keV
-        self.calibration_factor_b1 = 0.652
-        self.calibration_factor_b2 = 0.668
+        self.silizium = 3.6e-3
         self.my_path2irena = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/data/'
+        # mapping 
+        self.upper_trigger_real = {
+                'TO1' : 'TO2',
+                'TO2' : 'TO5',
+                'TO3' : 'TO9',
+                'TO4' : 'TO13',
+                'TO5' : 'TO16',
+                'TO6' : 'TO11',
+                'TO7' : 'TO7',
+                'TO8' : 'TO12'
+                }
+
+        self.lower_trigger = {
+                'TU2': 10,
+                'TU5': 13,
+                'TU9': 8,
+                'TU13': 7,
+                'TU16': 4,
+                'TU11': 14,
+                'TU7': 1, 
+                'TU12': 2
+                }
+        self.lower_trigger_real = {
+            'TU1' : 'TU2',
+            'TU2' : 'TU5',
+            'TU3' : 'TU9',
+            'TU4' : 'TU13',
+            'TU5' : 'TU16',
+            'TU6' : 'TU11',
+            'TU7' : 'TU7',
+            'TU8' : 'TU12'
+        }
+        self.bgo = {'B1':'B2','B2':'B1'}
+
         if delete:
             self.delete_old()
             self.download()
@@ -75,12 +113,21 @@ class Hist:
 
     def average(self):
         for column in self.data.columns[1:]:
-            kernel = np.ones(1)/1
+            kernel = np.ones(4)/4
             self.data[column] = np.convolve(self.data[column], kernel, mode='same')
 
     def plot_hist(self):
+        self.calibration_factors = pd.read_csv('fitting_values_gesamt.csv', delimiter=',')
         self.list_names_TO = ['TO1', 'TO2', 'TO3', 'TO4', 'TO5', 'TO6', 'TO7', 'TO8' ]
         self.list_names_TU = ['TU1', 'TU2', 'TU3', 'TU4', 'TU5', 'TU6', 'TU7', 'TU8' ]
+        colors = [
+            'blue', 'orange', 'green', 'red', 'purple',
+            'brown', 'pink', 'gray', 'olive', 'cyan',
+            'navy', 'darkorange', 'darkgreen', 'darkred', 'indigo',
+            'saddlebrown', 'lightcoral', 'dimgrey', 'darkolivegreen', 'teal'
+        ]
+   
+        color_iterator = cycle(colors)
         print('Sum of counts:', np.sum([np.sum(self.data[column]) for column in self.data.columns[1:]]))
         fig, ax = plt.subplots(figsize=(20,10))
         for column in self.data.columns[1:]:
@@ -101,29 +148,28 @@ class Hist:
                 #x,y = self.fit_gauss(self.data['mV'].to_numpy(),self.data[column].to_numpy(), init_vals=[10000000,150,40 ,   0.87,12.0,13283.0,   0.6, 0.2 ,714, 4000])
                 #ax.plot(x,y,linestyle='--',color='black')
                 pass
-            if column == 'B1' or column == 'B2' and column != 'B1' and column != 'B2':
-                #ax.plot(self.data['mV'], self.data[column], label=column)
-                #ax.legend()
-                pass
-            if column == 'B1':
-                ax.plot(self.data['mV']/ self.calibration_factor_b1, self.data[column], label=column)
-                print('B1 number of counts:', np.sum(self.data[column][self.data['mV']/ self.calibration_factor_b1 > 50]))
-            if column == 'B2':
-                ax.plot(self.data['mV']/ self.calibration_factor_b2, self.data[column], label=column)
-                print('B2 number of counts:', np.sum(self.data[column][self.data['mV']/ self.calibration_factor_b2 > 50]))
-                
+            if 'O' in column:
+                combined_dict = {**self.upper_trigger_real, **self.lower_trigger_real, **self.bgo}
+                ax.plot(self.data['mV']/self.silizium/ self.calibration_factors.loc[self.calibration_factors['name'] == column, 'e'].values[0], self.data[column], label = combined_dict[column], linewidth=2)
+            if 'U' in column:
+                combined_dict = {**self.upper_trigger_real, **self.lower_trigger_real, **self.bgo}
+                ax.plot(self.data['mV']/self.silizium/ self.calibration_factors.loc[self.calibration_factors['name'] == column, 'e'].values[0], self.data[column], label = combined_dict[column], linewidth=2, linestyle='--')
+            if 'B' in column:
+                combined_dict = {**self.upper_trigger_real, **self.lower_trigger_real, **self.bgo}
+                ax.step(self.data['mV']/self.silizium/ self.calibration_factors.loc[self.calibration_factors['name'] == column, 'e'].values[0], self.data[column], label = combined_dict[column], linewidth=2)
         # die Dioden sind unterschiedlich aufgeklebt
-        plt.title(f'Histogramm of {self.dat_name[:-4]}')
+        plt.title(f'Histogram of {self.dat_name[:-4]}')
         ax.set_yscale('log')
-        #ax.set_xlim(0,100)
-        ax.set_xlabel('mV')
+        ax.set_xlim(15000,65000)
+        ax.set_ylim(120,30000) 
+        ax.set_xlabel('emitted photoelectrons')
         ax.set_ylabel('Counts')
         ax.grid()
-        ax.legend()
+        ax.legend(loc='upper right', ncol = 2)
         # create folder for plots if it not exists already
         if not os.path.exists(f'Plots_{self.dat_name[:-4]}'):
             os.makedirs(f'Plots_{self.dat_name[:-4]}')
-        plt.savefig(f'Plots_{self.dat_name[:-4]}/' + f'B1_B2.pdf')
+        plt.savefig(f'Calibrated_histogramm_all.pdf')
         plt.show()
 
     def fit_landau(self, x_data, y_data):

plot_mpvs.py

@@ -6,9 +6,12 @@ import os
 # add picture of mapping
 import matplotlib.image as image
 from matplotlib.offsetbox import (OffsetImage, AnnotationBbox)
+# inrease font size
+import matplotlib as mpl
+mpl.rcParams['font.size'] = 15
 class mpv_s:
     def __init__(self):
-        self.target_directory = 'Results2023-12-21-AHBGOC-14-langzeit/mpvs.txt'
+        self.target_directory = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/mpvs.txt'
         # root of coordinate i mid of BGO crystal
         self.coordinates = pd.DataFrame({
             'names' : ['TO2_TU2', 'TO5_TU5', 'TO9_TU9', 'TO13_TU13', 'TO16_TU16', 'TO11_TU11', 'TO7_TU7', 'TO12_TU12'],
@@ -17,11 +20,13 @@ class mpv_s:
         })
         self.y_axis = ['TO2_TU2', 'TO5_TU5', 'TO9_TU9', 'TO13_TU13', 'TO16_TU16']
         self.x_axis = ['TO11_TU11', 'TO7_TU7', 'TO12_TU12']
+        self.silizium = 3.6e-3
         self.plot_xs()
         #self.plot_ys()
     def plot_xs(self):
         data = pd.read_csv(self.target_directory, delimiter=',')
-        fig, ax = plt.subplots(3,1,figsize=(5, 10))
+        fig, ax = plt.subplots(3,1,figsize=(15, 10))
+        position_diodes = 45
         for i in range(len(data['name'])):
             numbers = re.findall(r'\d+', data['name'][i])
             name = data['name'][i][:-5]
@@ -29,30 +34,40 @@ class mpv_s:
             if name in self.y_axis:
                 y_coordinate = self.coordinates['y'][self.coordinates['names'] == name]
                 # apply the same plots but with subplots
-                ax[0].errorbar(y_coordinate, data['b1b2_mpvs'][i], yerr=data['b1b2_errors'][i], fmt='o', label= f'coinc of TO{numbers[0]} and TU{numbers[1]}')
+                ax[0].errorbar(y_coordinate, data['b1b2_mpvs'][i], yerr=data['b1b2_errors'][i], fmt='o', label= f'coinc. of TO{numbers[0]} and TU{numbers[1]}')
                 # also b1 and b2
-                ax[1].errorbar(y_coordinate, data['b1_mpvs'][i], yerr=data['b1_errors'][i], fmt='o', label= f'coinc of TO{numbers[0]} and TU{numbers[1]}')
-                ax[2].errorbar(y_coordinate, data['b2_mpvs'][i], yerr=data['b2_errors'][i], fmt='o', label= f'coinc of TO{numbers[0]} and TU{numbers[1]}')
+                ax[1].errorbar(y_coordinate, data['b1_mpvs'][i], yerr=data['b1_errors'][i], fmt='o', label= f'coinc. of TO{numbers[0]} and TU{numbers[1]}')
+                ax[2].errorbar(y_coordinate, data['b2_mpvs'][i], yerr=data['b2_errors'][i], fmt='o', label= f'coinc. of TO{numbers[0]} and TU{numbers[1]}')
 
                 #legend,etc
-                ax[0].set_title("mpv's of seen energy by B1 and B2")
+                ax[0].set_title("mpvs of seen energy by B1 and B2")
                 ax[0].set_xlabel('y position in mm')
-                ax[0].set_ylabel('mpv[$keV_{Si}$]')
-                ax[0].legend()
+                ax[0].set_ylabel('mpv')
                 ax[0].grid()
 
-                ax[1].set_title("mpv's of seen energy by B1")
+
+                ax[1].set_title("mpvs of seen energy by B1")
                 ax[1].set_xlabel('y position in mm')
-                ax[1].set_ylabel('mpv[$keV_{Si}$]')
-                ax[1].legend()
+                ax[1].set_ylabel('mpv')
                 ax[1].grid()
 
-                ax[2].set_title("mpv's of seen energy by B2")
+                ax[2].set_title("mpvs of seen energy by B2")
                 ax[2].set_xlabel('y position in mm')
-                ax[2].set_ylabel('mpv[$keV_{Si}$]')
-                ax[2].legend()
+                ax[2].set_ylabel('mpv')
                 ax[2].grid()
         
+        # plot 2 vertical lines which represent the positions of the 2 readout diodes, also name that vertical line
+        ax[0].axvline(x=position_diodes, color='black', linestyle='--', linewidth=2, label='B1')
+        ax[0].axvline(x=-position_diodes, color='blue', linestyle='--', linewidth=2, label='B2')
+        ax[1].axvline(x=position_diodes, color='black', linestyle='--', linewidth=2, label='B1')
+        ax[1].axvline(x=-position_diodes, color='blue', linestyle='--', linewidth=2, label='B2')
+        ax[2].axvline(x=position_diodes, color='black', linestyle='--', linewidth=2, label='B1')
+        ax[2].axvline(x=-position_diodes, color='blue', linestyle='--', linewidth=2, label='B2')
+        # write B1 and B2 beside the vertical lines
+
+        ax[0].legend(ncol=2)
+        ax[1].legend(ncol=2)
+        ax[2].legend(ncol=2)
         # more space bewteen figures
         fig.tight_layout()
 
@@ -64,7 +79,7 @@ class mpv_s:
         ab = AnnotationBbox(imagebox, (0.2, 0.85), xycoords='figure fraction', frameon=False)
         fig.add_artist(ab)
         '''
-        plt.savefig('Results2023-12-21-AHBGOC-14-langzeit/mpvs.png')
+        plt.savefig('/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/mpvs.pdf', bbox_inches='tight')
         plt.show()
 
     def plot_ys(self):

round_DIN1333.py

@@ -0,0 +1,168 @@
+from decimal import Decimal as dec
+
+#Rundet nach DIN1333
+def round_DIN1333(value,error):
+    val_str = str(value)
+    err_str = str(error)
+
+# Fehlerrundung    
+    #Untersuche ob Zahl dezimal oder Zehnerpotenz und wandel in dezimal um
+    if err_str.find('e-') != -1:
+        c = err_str.find('e-')
+        err_str_dummy = '0.' + (int(err_str[c+2:])-1)*'0'
+        if '.' in err_str:
+            err_str_dummy += err_str[0] + err_str[2:c]
+        else:
+            err_str_dummy += err_str[:c]
+        err_str = err_str_dummy
+    elif err_str.find('e+') != -1:
+        c = err_str.find('e-')
+        if '.' in err_str:
+            err_str_dummy = err_str[0] + err_str[2:c]
+            err_str_dummy += (int(err_str[c+2:]) - len(err_str[2:c])) * '0'
+        else:
+            err_str_dummy = err_str[0] + int(err_str[c+2:]) * '0'
+        err_str = err_str_dummy
+            
+    
+    round_err = ''
+    # Finde Position des Punktes und merke Anzahl der Stellen vor und nach dem Punkt
+    point_i = err_str.find('.')
+    if point_i == -1:
+        bf_point = len(err_str)  
+    else:
+        af_point = len(err_str) - point_i - 1
+        bf_point = point_i
+    
+    #Finde erste von 0 verschiedene Stelle
+    round_i = 0
+    while round_i < len(err_str) and round_i != -1:
+        if round_i == len(err_str)-1 and err_str[round_i] == '0':
+            round_i = -1
+        elif err_str[round_i] not in ['.','0']:
+            break
+        else:
+            round_i += 1
+    
+    # Wenn die erste ungleich 0 Stelle kleiner der Position des Punktes ist, muss zwangsläufig auf die erste Stelle gerundet werden. 
+    # Außer diese ist eine 1 oder 2.
+    # Sonst wird aufgerundet        
+    if round_i <= abs(point_i):
+        if err_str[0] not in ['1','2']:
+            resolution=(1,0)
+            round_err = str(int(err_str[0])+1) + (bf_point-1)*'0'
+        elif err_str[0] in ['1','2']:
+            if bf_point != 1:
+                resolution = (2,0)
+                round_err = err_str[0] + str(int(err_str[1])+1) + (bf_point-2)*'0'
+            else:
+                resolution = (0,1)
+                round_err = err_str[:2] + str(int(err_str[2])+1)
+    
+    # Sonst wird nach dem Komma gerundet. Auch hier gilt bei 1 oder 2 wird auf die nächste Stelle gerundet.
+    # Wenn nicht 1 oder 2 wird die zurundene Stelle um 1 erhöht, da der Fehler immer aufgerundet wird.
+    else:
+        if err_str[round_i] not in ['1','2']:
+            resolution = (0,round_i-2)
+            round_err = str(dec(err_str[:round_i+1]) + dec('0.' + (round_i-2)*'0' + '1'))
+        else:
+            resolution = (0,round_i-1)
+            round_err = str(dec(err_str[:round_i+2]) + dec('0.' + (round_i-1)*'0' + '1'))
+    
+# Messwertrundung
+    round_val = ''
+    
+    bf_point_err = 0
+    af_point_err = 0
+    
+    point_err = round_err.find('.')
+
+
+    if point_err == -1:
+        bf_point_err = len(round_err)
+    else:
+        bf_point_err = point_err
+        af_point_err = len(round_err) - point_err - 1
+        
+        
+        
+    bf_point_val = 0
+    af_point_val = 0
+    
+    point_val = val_str.find('.')
+    if point_val == -1:
+        bf_point_val = len(val_str)
+    else:
+        bf_point_val = point_val
+        af_point_val = len(val_str) - point_val - 1
+           
+        
+    if af_point_err != 0 and af_point_err == af_point_val:
+        round_val = val_str
+    elif round_err[0] != 0:
+        #print('Case 1')
+        if point_err == -1:
+            i = round_err.find('0')
+            if i == -1:
+                index = bf_point_val - bf_point_err
+            else:
+                index = bf_point_val - bf_point_err + i -1
+            if bf_point_val > 1:
+                if val_str[index+1] != '.':
+                    if val_str[index+1] in ['0','1','2','3','4']:
+                        round_val = val_str[:index+1] + (len(round_err)-1)*'0'
+                    else:
+                        addition = int('1' + (len(round_err)-1)*'0')
+                        round_val = str(int(val_str[:index+1] + (len(round_err)-1)*'0') + addition)
+                else:
+                    if val_str[index+2] in ['0','1','2','3','4']:
+                        round_val = val_str[:index+1]
+                    else:
+                        round_val = str(int(val_str[:index+1]) + 1)
+            elif index == len(val_str)-3:
+                round_val = val_str[:bf_point_val]
+            else:
+                if bf_point_err - bf_point_val > 1:
+                    round_val = '0'
+                elif bf_point_err - bf_point_val == 1:
+                    if val_str[0] in ['0','1','2','3','4']:
+                        round_val = '0'
+                    else:
+                        round_val = '1' + (len(round_err)-1)*'0'
+        else:
+            if af_point_err >= af_point_val:
+                round_val = val_str + (len(round_err)-2 - af_point_val)*'0'
+            else:
+                round_val = val_str[:point_val+af_point_err]
+                index = point_val+af_point_err
+                if val_str[index+1] in ['0','1','2','3','4']:
+                    round_val += val_str[index]
+                elif val_str[index+1] in ['5','6','7','8','9']:
+                    addition = float('0.' + (af_point_err-1)*'0' + '1')
+                    #print(addition)
+                    round_val = str(float(val_str) + addition)[:index+1]
+                    
+    elif af_point_err > af_point_val:
+        round_val = val_str + (len(round_err)-2 - af_point_val)*'0'
+    elif af_point_err < af_point_val:
+        round_val = val_str[:point_val+af_point_err]
+        index = point_val+af_point_err
+        if val_str[index+1] in ['0','1','2','3','4']:
+            round_val += val_str[index]
+        elif val_str[index+1] in ['5','6','7','8','9']:
+            addition = float('0.' + (af_point_err-1)*'0' + '1')
+            #print(addition)
+            round_val = str(float(val_str) + addition)[:index+1]
+
+           
+    return round_val,round_err
+
+def multiply_with_error(self,value):
+     # Set precision to a sufficiently high value
+    getcontext().prec = 50  # You can adjust the precision as needed
+    # Convert values to Decimal
+    value_decimal = Decimal(str(value))
+    factor_decimal = Decimal(str(1000000))
+    # Perform multiplication with arbitrary precision
+    new_value = value_decimal * factor_decimal
+    return float(str(new_value).rstrip('0').rstrip('.'))

runde_fehler.py

@@ -0,0 +1,29 @@
+import math
+
+def round_propper(value, error):
+    value_decimal_position =find_decimal_position(value)
+    error_decimal_position =find_decimal_position(error)
+    shift = error_decimal_position - value_decimal_position
+    value = float_to_string_conversion(value)
+    error = float_to_string_conversion(error)
+    for position, number in enumerate(error):
+        if int(number) not in [0,1,2]:
+            pass
+        else:
+            if error[position + 1] != 0:
+                error[position] == str(int(error[position]) + 1)
+                round_error_position = position 
+                break
+    # now round the value at the same position, it must be relative to the comma
+    round_value_position = round_error_position - shift
+    if value[round_value_position +1] < 5:
+        pass
+    else:
+        value[round_value_position] == str(value[round_value_position] +1)
+def float_to_string_conversion(float_number):
+    return str(float_number).replace(".", "").lstrip("0")
+
+def find_decimal_position(float_number):
+    number_str = str(float_number)
+    return number_str.find(".")
+

test.py

@@ -0,0 +1,27 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+data = np.ones((4,4))
+x = np.array([1,2,3,4])
+y = np.array([1,2,3,4])
+# all rows from the columns 2 set to two, so the matrice displayed how it is shown here.
+data[:,2] = 2
+print(data)
+
+plt.pcolormesh(x,y,data)
+plt.show()
+
+data_x = [1,1,1,1,1,1]
+data_y = [1,1,1,0,0,0]
+
+hist, xedges, yedges = np.histogram2d(data_x, data_y, bins=2)
+print('-----')
+print(hist)
+print(xedges)
+print(yedges)
+plt.pcolormesh(xedges, yedges, hist)
+# colorbar
+plt.colorbar()
+plt.show()
+
+# x data stays x data, y data stays y data

test_fit.py

@@ -0,0 +1,35 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+import pandas as pd
+from scipy.optimize import curve_fit
+
+
+class test_fit:
+    def __init__(self) -> None:
+        self.file = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/1d_hist_muon/TO16_TU16m.txt'
+        self.fit_plot()
+
+    def fit_plot(self):
+        data = np.loadtxt(self.file, comments='#', delimiter=',')
+        bins = data[:,0]
+        histogramm = data[:,1]
+        mask = (bins > 50) & (bins < 400)
+        popt, pcov = curve_fit(self.mips, bins[mask], histogramm[mask], p0=[1, 60, 10])
+        print(popt)
+        histogramm_fit = self.mips(bins, *popt)
+        plt.yscale('log')
+        plt.xlim(0,200)
+        plt.ylim(1,100)
+        plt.plot(bins, histogramm, label='data')
+        plt.plot(bins, histogramm_fit, label='fit')
+        plt.show()
+
+    def mips(self,x, a,e,s):
+        # THE FUNCTION USED FOR FITTING
+        return a*self.landau((x-e)/s)
+    
+    def landau(self,l):
+        return np.sqrt(np.exp(-l-np.exp(-l))/2/np.pi)
+
+test_fit()