ADDED landau_sum.py AND landau_trigger.py

BigBGO.py

@@ -71,12 +71,12 @@ def main():
             data = [[],[],[],[],[],[],[]]
             for line in f:
                 l = line.split()
-                if len(l) > 22:
+                if not 'EOF' in l:
                     if l[0] == 'EI':
                         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 CmV(l,ch[3]) >= thr[3]):
                                 for i in range(0,len(ch)):
-                                    data[i].append(CmV(l,ch[i]))
+                                    data[i].append(CmV(l,ch[i])*u[i])
         h = do_hist(data)
         if args.plot:
             for i in range(0,len(ch)):
@@ -91,7 +91,7 @@ def main():
             data = [[],[]]
             for line in f:
                 l = line.split()
-                if len(l) > 22:
+                if not 'EOF' in l:
                     if l[0] == 'EI':
                         if not args.trigger or (args.trigger and CmV(l,ch[3]) >= thr[3]):
                             if not args.muon or (args.muon and ismuon(l)):  
@@ -132,11 +132,11 @@ def Cphase(line, i):
 def do_hist(data):
     bin = int((maxV-minV)/resV)
     hist, bins = np.histogram(data[0], bins=bin, range=(minV,maxV), density=False)
-    df = pd.DataFrame({'mV':bins[1:]/2+ bins[:-1]/2, name[0]: hist})
+    df = pd.DataFrame({'keV':bins[1:]/2+ bins[:-1]/2, name[0]: hist})
     for i in range(1, len(ch)):
         hist, bins = np.histogram(data[i], bins=bin, range=(minV,maxV), density=False)
-        df1 = pd.DataFrame({'mV':bins[1:]/2+ bins[:-1]/2, name[i]: hist})
-        df = pd.merge(df, df1 , left_on = 'mV', right_on = 'mV')
+        df1 = pd.DataFrame({'keV':bins[1:]/2+ bins[:-1]/2, name[i]: hist})
+        df = pd.merge(df, df1 , left_on = 'keV', right_on = 'keV')
     return df
 
 def do_hist_sum(data):

histplot.py

@@ -49,7 +49,7 @@ def main():
     
     channels = []
     end = '.hist'
-    keV = False
+    energy = 'mV'
     raw = args.raw
     u = args.norm
     if not args.channels == None:
@@ -62,7 +62,7 @@ def main():
     else:
         if not args.dhist == None:
             if 'keV' in args.dhist:
-                keV = True
+                energy = 'keV'
             for ch in args.dhist:
                 channels.append(str(ch))
             end = '.'+args.filepath.split('.')[-1]
@@ -84,9 +84,12 @@ def main():
         raw = True
     plt.figure(figsize=(17,12))
 
-    print_data(data,Itime,channels, u, keV)
+    if 'keV' in data:
+        energy = 'keV'
 
-    layout(tit,raw,args.limits,keV)
+    print_data(data,Itime,channels, u, energy)
+
+    layout(tit,raw,args.limits,energy)
 
 
     if args.export == None:
@@ -122,10 +125,7 @@ def read_data(file, end, raw, dhist):
 
     return data, Itime
 
-def print_data(data, Itime, channels, u,keV): 
-    energy = 'mV'
-    if keV:
-        energy = 'keV'
+def print_data(data, Itime, channels, u,energy): 
     resV = 0.838214/2 
     if not u == None:
         if len(u) == len(channels):
@@ -145,11 +145,11 @@ def print_data(data, Itime, channels, u,keV):
         
               
 
-def layout(tit,raw, limits,keV):
+def layout(tit,raw, limits,energy):
 
-    energy = 'mV'
-    if keV:
-        energy = 'keV'
+    # energy = 'mV'
+    # if keV:
+    #     energy = 'keV'
 
     fs = 15
     lim = [None, None, None, None]

landau_sum.py

@@ -0,0 +1,118 @@
+'''
+Autor: Ava Pohley
+'''
+
+import argparse
+import sys
+import matplotlib.pyplot as plt
+import pandas as pd
+import numpy as np
+from scipy import special
+from scipy import constants
+from scipy.optimize import curve_fit
+import os.path
+
+
+file = 'data/2023-07-24-BigBGO-Bi207h_sum'    
+keV = True
+
+tit = "Comparison sumA-sumB of '2023-04-07-BigBGO-Bi207'"
+plt.figure(figsize=(17,12))
+
+data = pd.read_csv(file+'.hist', sep = ' ', delim_whitespace=False, skiprows=0)
+    
+f = file
+Itime = None
+while len(f) > 0:
+    if os.path.isfile(f+'.Itime'):
+        Itime = int(pd.read_csv(f+'.Itime', sep = ' ', delim_whitespace=False, skiprows=0, names = ['time'])['time'][0])
+        f=''
+    else:
+        f=f[:-1]
+
+
+# FIT FUNCTIONS  
+def landau(x):
+    return np.sqrt(np.exp(np.subtract(-x,np.exp(-x)))/np.pi/2)
+
+def mips(x, a,e,s):   
+    return a * landau(np.subtract(x,e)/s)
+
+def find_par(data, energy, e1, e2):
+    a=0
+    go = True
+    while go == True:
+        a = a+1     
+        if data[energy][a]>=e1:
+            go = False 
+    go = True
+    b = a
+    while go == True:
+        b = b+1     
+        if data[energy][b]>=e2:
+            go = False      
+    return [a,b]   
+
+def find_par_new(f1, f2, pa,u, pos):
+    val_a = pa[1]-f1*pa[2]
+    val_b = pa[1]+f2*pa[2]
+    a = pos
+    b = pos
+    go = True
+    while go == True:
+        a = a-1     
+        if u[a]<=val_a:
+            go = False 
+    go = True
+    while go == True:
+        b = b+1     
+        if u[b]>=val_b:
+            go = False      
+    return [a,b]   
+
+a = 670
+b = 1000
+
+startparams = [3000,220,45]
+energy = 'keV'
+# plt.axvline(data[energy][a])
+# plt.axvline(data[energy][b])
+
+resV = 0.838214/2
+
+for ch in ['sumA','sumB']:   
+        sig = np.power(data[ch]+1, 0.5) 
+        par, cov = curve_fit(mips, data[energy][a:b], data[ch][a:b], startparams, sigma = sig[a:b], absolute_sigma=True)
+       
+        print(par)
+        perr = np.sqrt(np.diag(cov))
+        print(perr)
+    
+        if not Itime==None:
+            plt.plot(data[energy], np.multiply(data[ch], 3600/resV/Itime), label=ch+' ; $e$ = ('+ str(round(par[1],2)) + '$\pm$' +str(round(perr[1]+0.005,2))+')'+energy)
+            plt.plot(data[energy], np.multiply(mips(data[energy], *par), 3600/resV/Itime), 'k--') 
+        else:
+            plt.plot(data[energy], data[ch], label=ch+' ; $e$ = ('+ str(round(par[1],2)) + '$\pm$' +str(round(perr[1]+0.005,2))+')'+energy)
+            plt.plot(data[energy], mips(data[energy]*par), 'k--') 
+
+
+
+fs = 15
+ylab = 'counts [1/h/keV]'
+xlab = 'A [' + energy +']'
+plt.plot(0, 0,'k--', linewidth=1, label = 'landau_fit')
+plt.legend(fontsize = fs)
+plt.xticks(fontsize=fs)
+plt.yticks(fontsize=fs)
+plt.title(tit, fontsize = fs+10, pad = 15)
+
+plt.yscale('log')
+plt.ylabel(ylab, fontsize = fs+5)
+plt.xlabel(xlab, fontsize=fs+5)
+plt.ylim(3,100)
+plt.xlim(100,550) 
+
+plt.show()
+#plt.savefig('comparison_sums.pdf')
+
+#---------------------------------------------------------------------------------------------------------

landau_trigger.py

@@ -0,0 +1,120 @@
+'''
+Autor: Ava Pohley
+'''
+import matplotlib.pyplot as plt
+import pandas as pd
+import numpy as np
+from scipy import special
+from scipy import constants
+from scipy.optimize import curve_fit
+import os.path
+
+
+file = 'data/2023-07-24-BigBGO-Bi207h_muon_trigger'    
+keV = True
+
+tit = "Comparison single dioden when T1 has triggered '2023-04-07-BigBGO-Bi207'"
+plt.figure(figsize=(17,12))
+
+data = pd.read_csv(file+'.hist', sep = ' ', delim_whitespace=False, skiprows=0)
+ch = pd.read_csv(file+'.hist', sep = ' ', nrows=0)
+channels=['A1','A2','A3','B1','B2']
+    
+f = file
+Itime = None
+while len(f) > 0:
+    if os.path.isfile(f+'.Itime'):
+        Itime = int(pd.read_csv(f+'.Itime', sep = ' ', delim_whitespace=False, skiprows=0, names = ['time'])['time'][0])
+        f=''
+    else:
+        f=f[:-1]
+
+
+# FIT FUNCTIONS  
+def landau(x):
+    return np.sqrt(np.exp(np.subtract(-x,np.exp(-x)))/np.pi/2)
+
+def mips(x, a,e,s):   
+    return a * landau(np.subtract(x,e)/s)
+
+def find_par(data, energy, e1, e2):
+    a=0
+    go = True
+    while go == True:
+        a = a+1     
+        if data[energy][a]>=e1:
+            go = False 
+    go = True
+    b = a
+    while go == True:
+        b = b+1     
+        if data[energy][b]>=e2:
+            go = False      
+    return [a,b]   
+
+def find_par_new(f1, f2, pa,u, pos):
+    val_a = pa[1]-f1*pa[2]
+    val_b = pa[1]+f2*pa[2]
+    a = pos
+    b = pos
+    go = True
+    while go == True:
+        a = a-1     
+        if u[a]<=val_a:
+            go = False 
+    go = True
+    while go == True:
+        b = b+1     
+        if u[b]>=val_b:
+            go = False      
+    return [a,b]   
+
+a = [207,196,197,248,233]
+b = [270,270,270,320,320]
+energy = 'keV'
+startparams = [[150,85,10],[141,72,11],[141,72,11],[110,125,15],[110,110,16]]
+
+resV = 0.838214/2
+# plt.plot(data[energy], np.multiply(mips(data[energy], 150, 75, 15), 3600/resV/Itime),'--',label = '1')
+# plt.plot(data[energy], np.multiply(mips(data[energy], 140, 75, 15), 3600/resV/Itime),'--',label = '2')
+# plt.plot(data[energy], np.multiply(mips(data[energy], 130, 75, 15), 3600/resV/Itime),'--',label = '3')
+
+# plt.axvline(data[energy][a])
+# plt.axvline(data[energy][a1], color = 'r')
+
+for i in range (0,len(channels)):   
+        ch = channels[i]
+        sig = np.power(data[channels[i]]+1, 0.5) 
+        par, cov = curve_fit(mips, data[energy][a[i]:b[i]], data[ch][a[i]:b[i]], startparams[i], sigma = sig[a[i]:b[i]], absolute_sigma=True)
+       
+        print(par)
+        perr = np.sqrt(np.diag(cov))
+        print(perr)
+    
+        if not Itime==None:
+            #plt.plot(data[energy], np.multiply(data[channels[i]], 3600/resV/Itime), label=channels[i])
+            plt.plot(data[energy], np.multiply(data[ch], 3600/resV/Itime), label=ch+' ; $e$ = ('+ str(round(par[1],2)) + '$\pm$' +str(round(perr[1]+0.005,2))+')'+energy)
+            plt.plot(data[energy], np.multiply(mips(data[energy], *par), 3600/resV/Itime), 'k--') 
+        # else:
+        #     plt.plot(data[energy], data[ch], label=ch+' ; $e$ = ('+ str(round(par[1],2)) + '$\pm$' +str(round(perr[1]+0.005,2))+')'+energy)
+        #     plt.plot(data[energy], mips(data[energy]*par), 'k--') 
+
+
+
+fs = 15
+ylab = 'counts [1/h/keV]'
+xlab = 'A [' + energy +']'
+plt.plot(0, 0,'k--', linewidth=1, label = 'landau_fit')
+plt.legend(fontsize = fs)
+plt.xticks(fontsize=fs)
+plt.yticks(fontsize=fs)
+plt.title(tit, fontsize = fs+10, pad = 15)
+
+plt.yscale('log')
+plt.ylabel(ylab, fontsize = fs+5)
+plt.xlabel(xlab, fontsize=fs+5)
+plt.ylim(0.09,3)
+plt.xlim(0,300) 
+
+plt.show()
+#plt.savefig('comparison_sums.pdf')