Compare commits
No commits in common. "324a585e10d02ac3f117424e8753d43db6a1543f" and "19b87ce98fcc007da834bfa0033bb422f7851102" have entirely different histories.
324a585e10
...
19b87ce98f
24 changed files with 236 additions and 1453 deletions
BIN
DiodenBezeichnung.png
Normal file
BIN
DiodenBezeichnung.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 11 KiB |
|
|
@ -1,27 +0,0 @@
|
|||
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()
|
||||
|
|
@ -1,233 +0,0 @@
|
|||
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)
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -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[:-1]}',
|
||||
f'# measurement time in s: {self.measurement_time[:-2]}',
|
||||
f'# BGO1 energy deposit in mV: {i[3]}',
|
||||
f'# BGO2 energy deposit in mV: {i[4]}',
|
||||
f'# counts: {counts}\n'
|
||||
|
|
|
|||
|
|
@ -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 = 2e6
|
||||
self.chunksize = 3e6
|
||||
self.calibration_factor_b1 = 0.652
|
||||
self.calibration_factor_b2 = 0.668
|
||||
# to mV
|
||||
|
|
@ -65,7 +65,6 @@ 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
|
||||
|
|
|
|||
|
|
@ -1,197 +0,0 @@
|
|||
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)
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -51,7 +51,7 @@ neue Kanalzuordnung:
|
|||
|
||||
'''
|
||||
class Hist:
|
||||
def __init__(self, dat_name = '2023-12-11-AHBGOC-Bi207-08.dat' ):
|
||||
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/'
|
||||
|
|
@ -82,7 +82,6 @@ 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
|
||||
|
|
@ -141,16 +140,9 @@ 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
|
||||
if i == 'B1':
|
||||
else:
|
||||
# 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)
|
||||
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))
|
||||
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)
|
||||
|
|
@ -198,7 +190,7 @@ class Hist:
|
|||
|
||||
|
||||
source_file = '2023-12-21-AHBGOC-14-langzeit.dat'
|
||||
Hist()
|
||||
Hist(dat_name=source_file)
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,15 +0,0 @@
|
|||
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
|
||||
|
||||
|
|
@ -1,43 +0,0 @@
|
|||
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()
|
||||
|
|
@ -1,39 +0,0 @@
|
|||
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()
|
||||
|
|
@ -6,14 +6,11 @@ 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 = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/1d_hist_muon'
|
||||
self.path = '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],
|
||||
|
|
@ -42,7 +39,6 @@ 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 = []
|
||||
|
|
@ -75,65 +71,52 @@ 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
|
||||
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
|
||||
B2 = data[:, 1] / measurement_time
|
||||
B1 = data[:, 2] / measurement_time
|
||||
B1_B2 = data[:, 3] / measurement_time
|
||||
|
||||
bins_b1 = data[:, 0] / self.calibration_factor_b2
|
||||
bins_b2 = data[:, 0] / self.calibration_factor_b1
|
||||
bins_b1_b2 = data[:, 0]
|
||||
bins_b1_b2 = data[:, 0] / self.calibration_factor_b1
|
||||
# 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/self.silizium, B1, label=f'B1 mpv: m = '+ f'{mpv_1_round}$\pm${e_error_1_round}', where='mid')
|
||||
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')
|
||||
name = f'TO{i}_TU{i}_B2'
|
||||
ax.step(bins_b2/self.silizium, B2, label='B2 mpv: m = ' + f'{mpv_2_round}$\pm${e_error_2_round}', where='mid')
|
||||
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')
|
||||
name = f'TO{i}_TU{i}_B1+B2'
|
||||
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')
|
||||
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')
|
||||
|
||||
# Deactivate x-axis labels
|
||||
if i == self.right:
|
||||
ax.set_xticks([25000,50000,75000])
|
||||
ax.set_xticks([100,200,300])
|
||||
ax.grid(True)
|
||||
else:
|
||||
# Adjust x-axis ticks for compactness
|
||||
ax.set_xticks([25000,50000,75000, 100000])
|
||||
ax.set_xticks([100,200,300,400])
|
||||
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, top = 1)
|
||||
ax.set_xlim(0, 500/self.silizium)
|
||||
ax.set_ylim(bottom = 0.9/ measurement_time)
|
||||
ax.set_xlim(0, 500)
|
||||
ax.legend(loc='upper right', bbox_to_anchor=(1.2, 1.3))
|
||||
# Adding a text block
|
||||
text_block = 'x-axis: emitted photoelectrons\n y-axis: number of events in 1/h/mV'
|
||||
text_block = 'x-axis: Energydeposit diode per diode in $keV_{Si}$\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)
|
||||
|
||||
|
|
@ -150,8 +133,8 @@ class plot_1d_hist:
|
|||
plt.show()
|
||||
|
||||
def fit_landau(self, x_data, y_data, sum = False):
|
||||
mid_cut_low = 30
|
||||
mid_cut_high = 130
|
||||
mid_cut_low = 20
|
||||
mid_cut_high = 200
|
||||
if sum:
|
||||
mid_cut_low = 80
|
||||
mid_cut_high = 400
|
||||
|
|
@ -160,13 +143,12 @@ class plot_1d_hist:
|
|||
y_s = []
|
||||
mpv_s = []
|
||||
e_errors = []
|
||||
for lower_cut in range(0,20,2):
|
||||
for lower_cut in range(0,10,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)
|
||||
|
|
@ -180,13 +162,14 @@ 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 = x_data[:-1] + (x_data[1] - x_data[0]) / 2
|
||||
y_data = y_data[1:]
|
||||
# First iteration
|
||||
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))
|
||||
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]
|
||||
|
|
@ -204,7 +187,7 @@ class plot_1d_hist:
|
|||
return a*np.exp((e-x)/s)
|
||||
|
||||
|
||||
def lan(self,x, a,e,s):
|
||||
def lan(self,x, a,e,s,ab,sb):
|
||||
'''
|
||||
a = amplitude
|
||||
e = position
|
||||
|
|
@ -212,15 +195,6 @@ class plot_1d_hist:
|
|||
ab = amplitude bg
|
||||
sb = width bg
|
||||
'''
|
||||
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('.'))
|
||||
return self.mips(x, a,e,s) + self.bg(x, ab, e, sb)
|
||||
|
||||
plot_1d_hist()
|
||||
|
|
|
|||
|
|
@ -8,11 +8,10 @@ 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 = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/2dhistmuons'
|
||||
self.path = 'Results2023-12-21-AHBGOC-14-langzeit/2dhistmuons'
|
||||
self.file_names = self.file_names()
|
||||
self.mapping = pd.DataFrame({
|
||||
# vertikale linie von diode zu diode
|
||||
|
|
@ -39,7 +38,7 @@ class plot_1d_hist:
|
|||
})
|
||||
self.calibration_factor_b1 = 0.652
|
||||
self.calibration_factor_b2 = 0.668
|
||||
self.silizium = 3.6e-3
|
||||
self.energy_photon_keV = 3.6e-3
|
||||
self.used = [2,5,9,13,16,11,7,12]
|
||||
self.right = 7
|
||||
self.minV = -100
|
||||
|
|
@ -50,15 +49,6 @@ 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 = []
|
||||
|
|
@ -92,38 +82,39 @@ 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_b2
|
||||
self.bins_new_b2 = self.bins_new / self.calibration_factor_b1
|
||||
self.bins_new_b1 = self.bins_new / self.calibration_factor_b1
|
||||
self.bins_new_b2 = self.bins_new / self.calibration_factor_b2
|
||||
|
||||
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_b2/ self.silizium, self.bins_new_b1/ self.silizium, self.data, cmap=cmap, shading='auto', vmin=0.6)
|
||||
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/self.silizium, y/self.silizium, color='black',linewidth=2, linestyle='--',label='line with slope 1 and intercept 0')
|
||||
ax.plot(x, y, 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])
|
||||
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])
|
||||
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}')
|
||||
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}$')
|
||||
|
||||
plt.grid()
|
||||
plt.legend(loc = 'upper left')
|
||||
|
||||
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)
|
||||
|
||||
# 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.legend(loc='upper right', bbox_to_anchor=(1.1, 1.3))
|
||||
|
||||
# Adding a text block
|
||||
text_block = 'x-axis: B2 emitted photoelectrons\ny-axis: B1 emitted photoelectrons'
|
||||
text_block = 'x-axis: Energydeposit B2 in $keV_{Si}$\n y-axis: Energydeposit B1 in $keV_{Si}$'
|
||||
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)
|
||||
|
||||
|
|
@ -136,7 +127,8 @@ 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 = 400)
|
||||
plt.savefig('Overview_vertical_coinc_top_view_2d_hist.png', dpi = 300)
|
||||
plt.show()
|
||||
|
||||
def line_data_from_matrice(self):
|
||||
# i want to create histogram x and y data from the 2d histogram matrix
|
||||
|
|
@ -151,8 +143,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_b2[index_column])
|
||||
y.append(self.bins_new_b1[index_row])
|
||||
x.append(self.bins_new_b1[index_column])
|
||||
y.append(self.bins_new_b2[index_row])
|
||||
|
||||
self.x = np.array(x)
|
||||
self.y = np.array(y)
|
||||
|
|
|
|||
|
|
@ -1,39 +0,0 @@
|
|||
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()
|
||||
|
|
@ -8,15 +8,12 @@ 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': 25})
|
||||
plt.rcParams.update({'font.size': 20})
|
||||
'''
|
||||
This
|
||||
'''
|
||||
|
|
@ -28,7 +25,6 @@ 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 = []
|
||||
|
|
@ -36,121 +32,73 @@ 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
|
||||
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
|
||||
# corrected
|
||||
B2 = data[:, 2] / measurement_time
|
||||
B1 = data[:, 1] / measurement_time
|
||||
bins = data[:, 0]
|
||||
sum_B1_B2 = data[:, 3] / measurement_time / self.resV
|
||||
bins_b1 = bins / self.calibration_factor_b2
|
||||
bins_b2 = bins / self.calibration_factor_b1
|
||||
sum_B1_B2 = data[:, 3] / measurement_time
|
||||
bins_b1 = bins / self.calibration_factor_b1
|
||||
bins_b2 = bins / self.calibration_factor_b2
|
||||
numbers = re.findall(r'\d+', os.path.basename(file_name))
|
||||
# den threshold ermitteln
|
||||
if np.sum(B1) > 3:
|
||||
x_b1, y_b1, mpv_b1, e_error_b1, error, popt = self.fit_landau(bins_b1,B1)#
|
||||
|
||||
if np.sum(B1) > 10:
|
||||
x_b1, y_b1, mpv_b1, e_error_b1 = self.fit_landau(bins_b1,B1)
|
||||
name = fr'$m_{{\text{{TO{numbers[0]}TU{numbers[1]}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)
|
||||
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)
|
||||
|
||||
name = fr'$m_{{\text{{TO{numbers[0]}TU{numbers[1]}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)
|
||||
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)
|
||||
|
||||
name = fr'$m_{{\text{{TO{numbers[0]}TU{numbers[1]}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)
|
||||
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)
|
||||
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/ 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')
|
||||
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')
|
||||
else:
|
||||
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.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.set_yscale('log')
|
||||
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]))
|
||||
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]))
|
||||
plt.legend()
|
||||
plt.xlabel('emitted photoelectrons')
|
||||
plt.xlabel(r'Energy [$keV_{Si}$]')
|
||||
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.4)
|
||||
ab = AnnotationBbox(imagebox, (0.8, 0.45), xycoords='figure fraction', frameon=False)
|
||||
imagebox = OffsetImage(setup, zoom=0.3)
|
||||
ab = AnnotationBbox(imagebox, (0.8, 0.55), 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
|
||||
|
|
@ -164,9 +112,8 @@ class Plot1DHist:
|
|||
self.files = file_names_muon
|
||||
|
||||
def fit_landau(self, x_data, y_data, sum = False):
|
||||
mid_cut_low = 30
|
||||
# old cut 130
|
||||
mid_cut_high = 130
|
||||
mid_cut_low = 20
|
||||
mid_cut_high = 200
|
||||
if sum:
|
||||
mid_cut_low = 80
|
||||
mid_cut_high = 400
|
||||
|
|
@ -175,48 +122,39 @@ class Plot1DHist:
|
|||
y_s = []
|
||||
mpv_s = []
|
||||
e_errors = []
|
||||
error_list = []
|
||||
popt_list = []
|
||||
for lower_cut in range(0,20,2):
|
||||
for higher_cut in range(0,50,2):
|
||||
for lower_cut in range(0,10,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)
|
||||
#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)
|
||||
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(e_error)
|
||||
errors.append(error)
|
||||
e_errors.append(e_error)
|
||||
except Exception:
|
||||
pass
|
||||
# find best fit
|
||||
best_fit = np.argmin(errors)
|
||||
#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]
|
||||
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 = x_data[:-1] + (x_data[1] - x_data[0]) / 2
|
||||
y_data = y_data[1:]
|
||||
# First iteration
|
||||
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))
|
||||
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)
|
||||
# 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, error, popt
|
||||
return x_data, y, mpv, e_error
|
||||
|
||||
def mips(self,x, a,e,s):
|
||||
return a*self.landau((x-e)/s)
|
||||
|
|
@ -224,7 +162,11 @@ class Plot1DHist:
|
|||
def landau(self,l):
|
||||
return np.sqrt(np.exp(-l-np.exp(-l))/2/np.pi)
|
||||
|
||||
def lan(self,x, a,e,s):
|
||||
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
|
||||
|
|
@ -232,8 +174,7 @@ class Plot1DHist:
|
|||
ab = amplitude bg
|
||||
sb = width bg
|
||||
'''
|
||||
return self.mips(x, a,e,s)
|
||||
|
||||
return self.mips(x, a,e,s) + self.bg(x, ab, e, sb)
|
||||
|
||||
# Create an instance of the class
|
||||
source_file = '2023-12-21-AHBGOC-14-langzeit.dat'
|
||||
|
|
|
|||
|
|
@ -6,24 +6,14 @@ 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-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'
|
||||
def __init__(self, dat_name = '2023-12-06-AHBGOC-Bi207-07.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
|
||||
|
|
@ -33,11 +23,6 @@ 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)
|
||||
|
|
@ -46,44 +31,42 @@ 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,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)
|
||||
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))
|
||||
|
||||
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('fitting_values_6_B1.pdf')
|
||||
plt.savefig(file_path)
|
||||
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 'U232' in file_name or 'B' in file_name]
|
||||
file_names_muon = [file_name for file_name in file_names_muon if '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):
|
||||
|
|
@ -96,8 +79,8 @@ class Plot1DHist:
|
|||
pass
|
||||
|
||||
def fit_compton(self,x,y):
|
||||
mid_cut_low = 220
|
||||
mid_cut_high = 300
|
||||
mid_cut_low = 200
|
||||
mid_cut_high = 310
|
||||
init_vals = [3, 1, 260, 50, 100]
|
||||
|
||||
mask = (x > mid_cut_low) & (x < mid_cut_high)
|
||||
|
|
@ -107,17 +90,7 @@ 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,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('.'))
|
||||
return x,y, best_vals[2], error_e
|
||||
|
||||
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
|
||||
|
|
@ -132,6 +105,8 @@ 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))
|
||||
|
|
@ -140,7 +115,10 @@ 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))
|
||||
|
|
@ -202,4 +180,4 @@ class Plot1DHist:
|
|||
|
||||
# Create an instance of the class
|
||||
source_file = '2023-12-11-AHBGOC-Bi207-08.dat'
|
||||
plot_instance = Plot1DHist()
|
||||
plot_instance = Plot1DHist(dat_name = source_file)
|
||||
|
|
|
|||
|
|
@ -4,13 +4,11 @@ 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
|
||||
|
||||
|
|
@ -21,10 +19,8 @@ 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
|
||||
# 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.calibration_factor_b1 = 0.652
|
||||
self.calibration_factor_b2 = 0.668
|
||||
self.maxV = 3500
|
||||
self.resV = 0.838214/2
|
||||
self.bins = np.arange(self.minV, self.maxV, self.resV)
|
||||
|
|
@ -35,22 +31,9 @@ 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)
|
||||
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)
|
||||
|
|
@ -59,42 +42,32 @@ class Plot2DHist:
|
|||
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_b2
|
||||
self.bins_new_b2 = self.bins_new / self.calibration_factor_b1
|
||||
self.bins_new_b1 = self.bins_new / self.calibration_factor_b1
|
||||
self.bins_new_b2 = self.bins_new / self.calibration_factor_b2
|
||||
|
||||
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])]
|
||||
self.data = self.data[:, (self.bins[:-1] >= ylim[0]) & (self.bins[:-1] <= ylim[1])].T
|
||||
|
||||
# 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)
|
||||
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/ self.silizium, y/ self.silizium, color='black',linewidth=2, linestyle='--',label='line with slope 1 and intercept 0')
|
||||
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:
|
||||
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])
|
||||
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])
|
||||
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)
|
||||
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}$')
|
||||
|
||||
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('B2 emitted photoelectrons')
|
||||
plt.ylabel('B1 emitted photoelectrons')
|
||||
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 histogram coincidence of TO{numbers[0]} and TU{numbers[1]}')
|
||||
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
|
||||
|
|
@ -107,12 +80,10 @@ class Plot2DHist:
|
|||
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
|
||||
|
|
@ -122,13 +93,16 @@ 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.append(self.bins_new_b2[index_column])
|
||||
y.append(self.bins_new_b1[index_row])
|
||||
# x is B1 so
|
||||
x.append(self.bins_new_b1[index_column])
|
||||
y.append(self.bins_new_b2[index_row])
|
||||
|
||||
self.x = np.array(x)
|
||||
self.y = np.array(y)
|
||||
|
|
|
|||
|
|
@ -1,116 +0,0 @@
|
|||
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()
|
||||
|
|
@ -1,38 +0,0 @@
|
|||
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()
|
||||
86
plot_hist.py
86
plot_hist.py
|
|
@ -8,8 +8,6 @@ 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:
|
||||
|
|
@ -27,55 +25,19 @@ 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 = False, dat_name = '2023-12-21-AHBGOC-14-langzeit.dat'):
|
||||
def __init__(self, delete = True, 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.silizium = 3.6e-3
|
||||
self.calibration_factor_b1 = 0.652
|
||||
self.calibration_factor_b2 = 0.668
|
||||
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()
|
||||
|
|
@ -113,21 +75,12 @@ class Hist:
|
|||
|
||||
def average(self):
|
||||
for column in self.data.columns[1:]:
|
||||
kernel = np.ones(4)/4
|
||||
kernel = np.ones(1)/1
|
||||
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:]:
|
||||
|
|
@ -148,28 +101,29 @@ 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 '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)
|
||||
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]))
|
||||
|
||||
# die Dioden sind unterschiedlich aufgeklebt
|
||||
plt.title(f'Histogram of {self.dat_name[:-4]}')
|
||||
plt.title(f'Histogramm of {self.dat_name[:-4]}')
|
||||
ax.set_yscale('log')
|
||||
ax.set_xlim(15000,65000)
|
||||
ax.set_ylim(120,30000)
|
||||
ax.set_xlabel('emitted photoelectrons')
|
||||
#ax.set_xlim(0,100)
|
||||
ax.set_xlabel('mV')
|
||||
ax.set_ylabel('Counts')
|
||||
ax.grid()
|
||||
ax.legend(loc='upper right', ncol = 2)
|
||||
ax.legend()
|
||||
# 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'Calibrated_histogramm_all.pdf')
|
||||
plt.savefig(f'Plots_{self.dat_name[:-4]}/' + f'B1_B2.pdf')
|
||||
plt.show()
|
||||
|
||||
def fit_landau(self, x_data, y_data):
|
||||
|
|
|
|||
45
plot_mpvs.py
45
plot_mpvs.py
|
|
@ -6,12 +6,9 @@ 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 = '/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/mpvs.txt'
|
||||
self.target_directory = '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'],
|
||||
|
|
@ -20,13 +17,11 @@ 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=(15, 10))
|
||||
position_diodes = 45
|
||||
fig, ax = plt.subplots(3,1,figsize=(5, 10))
|
||||
for i in range(len(data['name'])):
|
||||
numbers = re.findall(r'\d+', data['name'][i])
|
||||
name = data['name'][i][:-5]
|
||||
|
|
@ -34,40 +29,30 @@ 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("mpvs of seen energy by B1 and B2")
|
||||
ax[0].set_title("mpv's of seen energy by B1 and B2")
|
||||
ax[0].set_xlabel('y position in mm')
|
||||
ax[0].set_ylabel('mpv')
|
||||
ax[0].set_ylabel('mpv[$keV_{Si}$]')
|
||||
ax[0].legend()
|
||||
ax[0].grid()
|
||||
|
||||
|
||||
ax[1].set_title("mpvs of seen energy by B1")
|
||||
ax[1].set_title("mpv's of seen energy by B1")
|
||||
ax[1].set_xlabel('y position in mm')
|
||||
ax[1].set_ylabel('mpv')
|
||||
ax[1].set_ylabel('mpv[$keV_{Si}$]')
|
||||
ax[1].legend()
|
||||
ax[1].grid()
|
||||
|
||||
ax[2].set_title("mpvs of seen energy by B2")
|
||||
ax[2].set_title("mpv's of seen energy by B2")
|
||||
ax[2].set_xlabel('y position in mm')
|
||||
ax[2].set_ylabel('mpv')
|
||||
ax[2].set_ylabel('mpv[$keV_{Si}$]')
|
||||
ax[2].legend()
|
||||
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()
|
||||
|
||||
|
|
@ -79,7 +64,7 @@ class mpv_s:
|
|||
ab = AnnotationBbox(imagebox, (0.2, 0.85), xycoords='figure fraction', frameon=False)
|
||||
fig.add_artist(ab)
|
||||
'''
|
||||
plt.savefig('/home/tomruge/Schreibtisch/UNI/Bachelorarbeit/PlotsAndData/Results2023-12-21-AHBGOC-14-langzeit/mpvs.pdf', bbox_inches='tight')
|
||||
plt.savefig('Results2023-12-21-AHBGOC-14-langzeit/mpvs.png')
|
||||
plt.show()
|
||||
|
||||
def plot_ys(self):
|
||||
|
|
|
|||
168
round_DIN1333.py
168
round_DIN1333.py
|
|
@ -1,168 +0,0 @@
|
|||
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('.'))
|
||||
|
|
@ -1,29 +0,0 @@
|
|||
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(".")
|
||||
|
||||
27
test.py
27
test.py
|
|
@ -1,27 +0,0 @@
|
|||
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
|
||||
35
test_fit.py
35
test_fit.py
|
|
@ -1,35 +0,0 @@
|
|||
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()
|
||||
Loading…
Add table
Add a link
Reference in a new issue