AHBGOx.py

@@ -1,218 +0,0 @@
-import numpy as np
-import argparse
-import fileinput
-import sys
-import pandas as pd
-from datetime import datetime
-import matplotlib.pyplot as plt
-
-parser = argparse.ArgumentParser(description="This script parses EI-data for CHAOS_jr and return histogramms")
-parser.add_argument("-v", "--verbose", choices = ['0','1','2'], help = "when zero only emit serious error messages; when one also emit warnings, when two emit informational messages")
-parser.add_argument("file", type = str, help="path/file were data is stored in")
-parser.add_argument('-hist', action='store_true', help = "makes histogramm of all channels")  # on/off flag
-parser.add_argument("-sum", action='store_true', help = "makes histogramms of sum(s) of BGOdiodes")  # on/off flag
-parser.add_argument('-Itime', action='store_true', help = "if chosen, Itime is calculatet and printet to filename.Itime")  # on/off flag
-parser.add_argument("-mV", action='store_true', help = "if chosen, data will be kept in mV changed to keV otherwise")  # on/off flag
-parser.add_argument('-muon', action='store_true', help = "only select muon data")  # on/off flag
-parser.add_argument("-xray", action='store_true', help = "only select xray data")  # on/off flag
-parser.add_argument("-trigger", type=str, nargs='+', help = "letters of Triggerdiodes to chose")
-parser.add_argument("-two", action='store_true', help = "if chosen, functions function for two diodes A1 and A2 if X=AHBGOD")  # on/off flag
-parser.add_argument("-zero", action='store_true', help = "if chosen, data will be normalized to 0°C, chose only for sum")  # on/off flag
-parser.add_argument("-export", type=str, nargs='+', help = "filpath to store output if other then currend directory/data/filename")
-args = parser.parse_args()
-
-file = args.file
-filename = file.split(".")[0]
-
-if args.verbose == None:
-    verbosity = 1
-else:
-    verbosity = args.verbose
-
-mV = 14000
-minV = -100
-maxV = 3500
-resV = 0.838214
-if args.trigger==None:
-    resV=resV/2
-pol = [1.65857396e-04, -7.56891495e-03, -3.66207616e+00,  3.70051587e+02] #polynom zur Temperaturumrechung
-
-global x
-x = None
-if 'AHBGOD' in file:
-    x = 'D'
-elif 'AHBGOC' in file:
-    x = 'C'
-
-print(x)
-
-def main():
-    if args.Itime:
-        Itime = find_Itime()
-        fname = filename + '.Itime'
-        outfile = open(fname, 'w')
-        outfile.write(str(Itime))
-        outfile.close()
-
-    else:
-        global thr
-        global ch
-        ch,name,thr,u = channels()
-
-        if name[0] == None and verbosity>0:
-            print("no definition selected. Chose -x C or -x D")
-            sys.exit()
-
-        global trig
-        trig=[]
-        nameadd='x'
-        if args.sum: nameadd = nameadd + '_sum'
-        if args.xray: nameadd = nameadd + '_xray'
-        if args.muon: nameadd = nameadd + '_muon'
-        if args.trigger:
-            nameadd = nameadd + '_trigger'
-            for t in args.trigger:
-                trig.append(t) 
-                nameadd = nameadd + '_' + str(t)
-        if args.zero: nameadd = nameadd + '_zero'
-        if args.mV: 
-            nameadd = nameadd + '_mV'
-            u = [1.0]*18
-
-        global trigchan 
-        trigchan = []
-        for t in trig:
-            for i in range(len(name)):
-                if t == name[i]:
-                    trigchan.append(i)
-        print(trigchan)
-
-        bin = int((maxV-minV)/resV)
-        hist = np.zeros((bin,18+1))
-        for i in range(0,bin):
-            hist[i][0]=minV + resV*i
-
-        with fileinput.input(files=(file), encoding="utf-8") as f: 
-            T='unknown'
-            for line in f:
-                l = line.split()
-                if not 'EOF' in l:
-                    if l[0] == 'H':
-                        T = adc_to_T(int(l[8]))
-                    elif l[0] == 'EI' and not T=='unknown':    
-                        if conditions(l,args):  
-                            for i in range(18):
-                                x = CmV(l,ch[i])*u[i]   
-                                if args.zero:
-                                    x = x*polynom(0,*pol)/polynom(T,*pol)
-                                if x>=minV and x<=maxV:
-                                    xx = int(x-minV/resV)
-                                    hist[xx,i+1] += 1
-
-        frame = pd.DataFrame(hist, columns = ['mV'] + name)
-        frame.to_csv(filename +nameadd+'.hist', sep = ' ', index = False)
-
-def channels():
-    ch = np.arange(18)
-    name=[None]*18
-    thr=[12]*18
-    u = [1.0]*18
-    if x == 'D':
-        name[0] = "A1"; ch[0] = 8; u[0] = 1/0.689
-        name[1] = "A2"; ch[1] = 2; u[1] = 1/0.680
-        name[2] = "A3"; ch[2] = 4; u[2] = 1/0.692
-        name[3] = "T1"; ch[3] = 7; u[3] = 1/0.68 #unsicher
-        name[4] = "B1"; ch[4] = 16; u[4] = 1/0.674
-        name[5] = "B2"; ch[5] = 13; u[5] = 1/0.668
-
-    elif x == 'C':
-        name[0] = "B1"; ch[0] = 17; thr[0] = 8; u[0] = 1/0.7
-        name[1] = "B2"; ch[1] = 16; thr[1] = 8; u[1] = 1/0.7
-        name[2] = "TO1"; ch[2] = 0; u[2] = 1/0.7
-        name[3] = "TO2"; ch[3] = 3; u[3] = 1/0.7
-        name[4] = "TO3"; ch[4] = 5; u[4] = 1/0.7
-        name[5] = "TO4"; ch[5] = 6; u[5] = 1/0.7
-        name[6] = "TO5"; ch[6] = 9; u[6] = 1/0.7
-        name[7] = "TO6"; ch[7] = 11; u[7] = 1/0.7
-        name[8] = "TO7"; ch[8] = 12; u[8] = 1/0.7
-        name[9] = "TO8"; ch[9] = 15; u[9] = 1/0.7
-        name[10] = "TU1"; ch[10] = 1; u[10] = 1/0.7
-        name[11] = "TU2"; ch[11] = 2; u[11] = 1/0.7
-        name[12] = "TU3"; ch[12] = 4; u[12] = 1/0.7
-        name[13] = "TU4"; ch[13] = 7; u[13] = 1/0.7
-        name[14] = "TU5"; ch[14] = 8; u[14] = 1/0.7
-        name[15] = "TU6"; ch[15] = 13; u[15] = 1/0.7
-        name[16] = "TU7"; ch[16] = 10; u[16] = 1/0.7
-        name[17] = "TU8"; ch[17] = 14; u[17] = 1/0.7
-
-    return ch, name, thr, u 
-
-def find_Itime():
-    Itime = 0
-    lasttime = 0
-    with fileinput.input(files=(file), encoding="utf-8") as f:
-        for line in f:
-            if line[0] == 'H':
-                time = int(line.split()[1])
-                diff = time - lasttime
-                if diff>0 and diff <= 120:
-                    Itime += diff
-                lasttime=time
-    return Itime
-
-def CmV(line, i):
-    return int(line[3*i+5])/mV
-
-def polynom(x,  a,b,c,d):
-    return np.multiply(np.multiply(np.multiply(a,x),x),x)+np.multiply(np.multiply(b,x),x)+np.multiply(c,x)+d
-
-def adc_to_T(adc):
-    R25 = 10000
-    R1 = 10000
-    b = 3940
-    dac = 62690
-    T = 0
-    alpha = np.log(R1/R25 *(dac/(dac-adc)-1))
-    T = 1/(alpha/b+1/298.15) -273.15
-    return T 
-
-
-def conditions(l,args):
-    if (not args.xray and not args.muon) or (args.xray and isXray(l)) or (args.muon and isMuon(l)):
-        if not args.trigger or (args.trigger and isTrigger(l)):
-            return True
-    else: return False
-
-def isXray(l):
-    if x == 'D':
-        if (CmV(l,ch[0]) >= thr[0] and CmV(l,ch[1]) < thr[1] and CmV(l,ch[2]) < thr[2]) or \
-        (CmV(l,ch[1]) >= thr[1] and CmV(l,ch[0]) < thr[0] and CmV(l,ch[2]) < thr[2]) or \
-        (CmV(l,ch[2]) >= thr[2] and CmV(l,ch[1]) < thr[1] and CmV(l,ch[0]) < thr[0]) or \
-        (CmV(l,ch[4]) >= thr[5] and CmV(l,ch[5]) < thr[5]) or \
-        (CmV(l,ch[5]) >= thr[6] and CmV(l,ch[4]) < thr[4]) or \
-        (CmV(l,ch[3]) >= thr[3] and CmV(l,ch[4]) < thr[4] and CmV(l,ch[5]) < thr[5]): return True
-    elif x == 'C':
-        if CmV(l,ch[0]) <= thr[0] or CmV(l,ch[1]) <= thr[1]: return True
-    else: return False
-
-def isMuon(l):
-    if x == 'D':
-        if not args.two and CmV(l,ch[4]) >= thr[4] and CmV(l,ch[5]) >= thr[5] and \
-        CmV(l,ch[0]) >= thr[0] and CmV(l,ch[1]) >= thr[1] and CmV(l,ch[2]) >= thr[2]: return True
-        if args.two and CmV(l,ch[4]) >= thr[4] and CmV(l,ch[5]) >= thr[5] and \
-        (CmV(l,ch[1]) >= thr[1] and CmV(l,ch[2]) >= thr[2]): return True
-        #CmV(l,A1) >= thr[0] and CmV(l,A2) >= thr[1]: return True
-    elif x == 'C':
-        if CmV(l,17) >= thr[0] and CmV(l,16) >= thr[1]: return True
-    else: return False
-
-def isTrigger(l):
-    if x == 'D':
-        if CmV(l,ch[3]) >= thr[3]: return True
-    elif x == 'C':
-        for t in trigchan:
-            if CmV(l,ch[t]) <= thr[t]: return False
-        return True
-    else: return False
-
-main()

BigBGO.py

@@ -96,7 +96,7 @@ def main():
                             if not args.trigger or (args.trigger and CmV(l,ch[3]) >= thr[3]):
                                 for i in range(0,len(ch)):
                                     if args.zero:
-                                        data[i].append(CmV(l,ch[i])*u[i]*polynom(0,*pol)/polynom(T,*pol))
+                                        data[i].append(CmV(l,ch[i])*u[i])*polynom(0,*pol)/polynom(T,*pol)
                                     else:
                                         data[i].append(CmV(l,ch[i])*u[i])
         h = do_hist(data)

landau_sum.py

@@ -16,11 +16,13 @@ from fit_functions import mips
 from fit_functions import find_par
 from fit_functions import find_par_new
 
+plt.figure(figsize=(15,10))
+
 file = 'data/2023-07-24-BigBGO-Bi207h_muon_sum'  
 tit = "Comparison sumA(3)-sumB(2) of '2023-07-24-BigBGO-Bi207'"  
 
-#file = 'data/2023-07-29-BigBGO-Bi207-2h_muon_sum'    
+file = 'data/2023-07-29-BigBGO-Bi207-2h_muon_sum'    
-#tit = "Comparison sumA(3)-sumB(2) of '2023-07-29-BigBGO-Bi207-2'"
+tit = "Comparison sumA(3)-sumB(2) of '2023-07-29-BigBGO-Bi207-2'"
 a = 610
 b = 1100
 
@@ -28,7 +30,7 @@ b = 1100
 # tit = "Comparison sumA(2 120°)-sumB(2) of '2023-08-07-BigBGO-3'"
 
 keV = True
-fig, ax = plt.subplots(figsize=(15,10))
+plt.figure(figsize=(17,12))
 
 data = pd.read_csv(file+'.hist', sep = ' ', delim_whitespace=False, skiprows=0)
     
@@ -63,33 +65,23 @@ for ch in ['sumA','sumB']:
             plt.plot(data[energy][a:b], lan(data[energy]*par)[a:b], 'k--') 
 
 
+
 fs = 15
-ax2=ax.twiny()
+ylab = 'counts [1/h/$keV_{Si}$]'
-ylab = 'counts [1/h/keV$_{Si}$]'
+xlab = 'A [' + energy +'$_{Si}$]'
-xlab = 'signal [' + energy +'$_{Si}$]'
+plt.plot(0, 0,'k--', linewidth=1, label = 'landau_fit')
-ax.plot(0, 0,'k--', linewidth=1, label = 'landau_fit')
+plt.legend(fontsize = fs)
-ax.legend(fontsize = fs)
+plt.xticks(fontsize=fs)
-ax.tick_params(axis='x', labelsize=fs)
+plt.yticks(fontsize=fs)
-ax2.tick_params(axis='x', labelsize=fs-5)
+plt.title(tit, fontsize = fs+10, pad = 15)
-ax.tick_params(axis='y', labelsize=fs)
-plt.title(tit, fontsize = fs+10, pad = 18)
 
 plt.yscale('log')
-ax.set_ylabel(ylab, fontsize = fs+5)
+plt.ylabel(ylab, fontsize = fs+5)
-ax.set_xlabel(xlab, fontsize=fs+5)
+plt.xlabel(xlab, fontsize=fs+5)
-#ax.set_ylim(0.3,100) #all
+#plt.ylim(5,100) #1
-#ax.set_xlim(50,850) #all
+plt.ylim(2,30) #2
-#ax.set_ylim(0.7,20) #shifted
+#plt.ylim(0.9,50) #3
-#ax.set_xlim(170,770) #shifted
+plt.xlim(100,450) 
-ax.set_ylim(0.9,30) #straight
-ax.set_xlim(100,550) #straight
-#ax.set_ylim(0.3,30) #120°
-#ax.set_xlim(100,800) #120°
-
-a, b = ax.get_xlim()
-
-ax2.set_xlim(a/0.00362,b/0.00362)
-ax2.set_xlabel('signal [counts of $e^{-}_{ph}$]', fontsize=fs)
 
 #plt.show()
 plt.savefig('plots/landau_'+file.split('/')[-1])

temperature_dependence.py

@@ -1,94 +0,0 @@
-import matplotlib.pyplot as plt
-import matplotlib.colors
-import numpy as np
-import math 
-import pandas as pd
-from datetime import date
-from datetime import datetime
-from scipy.optimize import curve_fit
-
-def polynom(x,  a,b,c,d):
-    return np.multiply(np.multiply(np.multiply(a,x),x),x)+np.multiply(np.multiply(b,x),x)+np.multiply(c,x)+d
-
-#fitparam plot_landau_for_temp_2_SUM(keV)
-a = [174.07, 187.61, 170.85, 154.53, 138.44, 126.15, 113.43, 104.51, 99.9, 97.58, 86.16, 85.08]
-err_a = [3.13, 2.93, 2.65, 2.4, 2.19, 2.01, 1.84, 1.69, 1.65, 1.6, 1.47, 1.47]
-e = [209.63, 239.4, 277.56, 315.6, 352.67, 384.33, 418.96, 443.88, 474.5, 494.25, 510.59, 516.59]
-err_e = [1.03, 1.14, 1.28, 1.4, 1.58, 1.74, 1.81, 2.1, 2.19, 2.15, 2.4, 2.5]
-s = [32.19, 38.66, 42.55, 47.57, 51.57, 56.56, 60.63, 65.99, 67.15, 70.15, 72.99, 73.12]
-err_s = [0.57, 0.59, 0.65, 0.73, 0.82, 0.91, 0.98, 1.09, 1.14, 1.16, 1.27, 1.31]
-
-a2 = [295.92, 174.89, 160.28, 150.87, 133.77, 119.64, 117.84, 107.74, 90.32, 95.87, 89.36, 69.86]
-err_a2 = [4.2, 2.88, 2.65, 2.42, 2.17, 1.98, 1.89, 1.77, 1.61, 1.59, 1.49, 1.37]
-e2 = [204.96, 231.3, 269.8, 307.33, 339.57, 377.43, 406.54, 439.53, 462.99, 485.35, 505.95, 520.52]
-err_e2 = [0.63, 1.17, 1.28, 1.43, 1.6, 1.74, 1.75, 1.98, 2.23, 2.26, 2.41, 2.64]
-s2 = [29.6, 37.44, 40.39, 45.97, 50.64, 54.81, 59.67, 62.57, 64.76, 69.47, 74.07, 70.68]
-err_s2 = [0.38, 0.61, 0.66, 0.74, 0.82, 0.91, 0.95, 1.04, 1.18, 1.19, 1.27, 1.44]
-
-TD = [44.43, 34.97, 25.46, 15.97, 6.68, -2.94, -12.35, -21.5, -30.48, -38.96, -46.87, -52.03]
-TD2 = [44.82, 35.45, 25.99, 16.48, 6.99, -2.6, -12.01, -21.22, -30.26, -38.9, -47.07, -51.54]
-
-fig, ax = plt.subplots(figsize=(17,12))
-
-par, cov = curve_fit(polynom, TD, e, sigma = err_e, absolute_sigma=True)
-par2, cov2 = curve_fit(polynom, TD2, e2, sigma = err_e2, absolute_sigma=True)
-
-perr = np.sqrt(np.diag(cov))
-perr2 = np.sqrt(np.diag(cov2))
-###############################################################
-n = polynom(0,*par)
-print(n)
-n2 = polynom(0,*par2)
-
-par3 = []
-perr3 = []
-for i in range(0,len(par)):
-    pa = np.array([par[i],par2[i]])
-    σpa = np.array([perr[i],perr2[i]])
-    g = 1/σpa**2
-    pa_mean = sum(g*pa)/sum(g)
-    σ2pa_mean = sum(g*(pa-pa_mean)**2)/(pa.shape[0]-1)/sum(g)
-    par3.append(pa_mean)
-    perr3.append(σ2pa_mean)
-
-
-print(par)
-print(perr)
-print(par2)
-print(perr2)
-print(par3)
-print(perr3)
-
-# Die eintreffenden teilchen deponieren JEDESMAL 18MeV im BGO, die dann das Licht erzeugen (Licht Temperaturabhängig)
-# 3.62 eV ist der Fano factor F, mean energy per electron-hole pair e-h in the high-energy limit, band-gap energy Eg the ratio A of phonon-carrier to carrier-carrier scattering, the optical phonon energies ~ωo, and the plasmon energy ~ωp l. 
-# keV / 0.00362 keV/e-   / 18MeV
-
-factor = 1/0.00362/18
-
-x = np.linspace(-52.0, 45.0, num=100)
-
-# #keV
-ax.errorbar(TD, np.multiply(e,factor), np.multiply(err_e,factor), elinewidth=2, capsize=5, fmt='bx', linewidth = 2, ecolor='b',markeredgecolor='b')
-ax.errorbar(TD2, np.multiply(e2,factor), np.multiply(err_e2,factor), elinewidth=2, capsize=5, fmt='rx', linewidth = 2, ecolor='r',markeredgecolor='r')
-ax.plot(TD, np.multiply(polynom(TD,*par),factor),'b--',linewidth = 1.8)
-ax.plot(TD2, np.multiply(polynom(TD2,*par2),factor),'r--',linewidth = 1.8)
-ax.plot(TD, np.multiply(polynom(TD,*par),factor),'b--',linewidth = 1.8,label = '$e(T)_{decreasing}$ = (' +  str(round(par[0]*10000*factor,1)) + '$ \cdot 10^{-4}$ $1/°C^3 $ $ T^3$ - ' + str(round(par[1]*-1000*factor,1)) + '$\cdot 10^{-3}$ $1/°C^2}$ $T^2$ - '+ str(round(par[2]*-1*factor,1)) + ' 1/°C $\cdot$ T + '+ str(round(par[3]*factor,1))+ ')  $e_{ph}/MeV$')
-ax.plot(TD2, np.multiply(polynom(TD2,*par2),factor),'r--',linewidth = 1.8,label = '$e(T)_{increasing}$  = (' + str(round(par2[0]*10000*factor,1)) + '$ \cdot 10^{-4}$ $1/°C^3 $ $T^3$ - ' + str(round(par2[1]*-1000*factor,1)) + '$\cdot 10^{-3}$ $1/°C^2$ $ T^2$ - '+ str(round(par2[2]*-1*factor,1)) + ' 1/°C $\cdot$ T + '+ str(round(par2[3]*factor,1))+ ')  $e_{ph}/MeV$')
-ax.plot(x, np.multiply(polynom(x,*par3),factor),'k-',linewidth = 1.8,label = '$e(T)_{mean}$       =  (' + str(round(par3[0]*10000*factor,1)) + '$ \cdot 10^{-4}$ $1/°C^3 $ $T^3$ - ' + str(round(par3[1]*-1000*factor,1)) + '$\cdot 10^{-3}$ $1/°C^2$ $ T^2$ - '+ str(round(par3[2]*-1*factor,1)) + ' 1/°C $\cdot$ T + '+ str(round(par3[3]*factor,1))+ ')  $e_{ph}/MeV$')
-
-print('$e(T)_{mean}$       =  (' + str(round(par3[0]*10000*factor,1)) + '$ \cdot 10^{-4}$ $1/°C^3 $ $T^3$ - ' + str(round(par3[1]*-1000*factor,1)) + '$\cdot 10^{-3}$ $1/°C^2$ $ T^2$ - '+ str(round(par3[2]*-1*factor,1)) + ' 1/°C $\cdot$ T + '+ str(round(par3[3]*factor,1))+ ')  $e_{ph}/MeV$')
-
-ax.legend(fontsize=15)
-ax.set_xlabel("T [°C]", fontsize=20)
-
-plt.title("Temperature dependence of light output in BGO", fontsize=24, pad =15)
-ax.set_ylabel("light output [$e_{ph}/MeV$]", fontsize=20)
-
-ax.tick_params(axis='x', labelsize=15)
-ax.tick_params(axis='y', labelsize=15)
-plt.yticks(fontsize=15)
-plt.ylim(3000,8600)
-#plt.grid(True)
-#plt.show()
-plt.savefig("temperature_dependence.pdf")
-