ADD CHAOShk.py ADD histplot2d.py

CHAOShk.py

@@ -0,0 +1,469 @@
+import numpy as np
+import argparse
+import fileinput
+import sys
+import pandas as pd
+from datetime import datetime
+from pathlib import Path
+
+parser = argparse.ArgumentParser(description="This script parses EI-data for AHBGO(C/D/S) and returns 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('-Itime', action='store_true', help = "if chosen, Itime is calculatet and printet to filename.Itime")  # on/off flag
+
+parser.add_argument('-hist', action='store_true', help = "makes histogramm of all channels using EI-lines")  # on/off flag
+parser.add_argument("-mV", action='store_true', help = "needs option -hist; if chosen, data will be kept in mV changed to keV otherwise")  # on/off flag
+parser.add_argument("-trigger", type=str, nargs='+', help = "needs option -EI; letters of Triggerchannels to chose")
+parser.add_argument("-xray", action='store_true', help = "needs option -EI; only select xray data")  # on/off flag
+
+parser.add_argument('-P', action='store_true', help = "evaluates pressure and temperature data from P-lines")  # on/off flag
+parser.add_argument('-C64', action='store_true', help = "evaluates counts per second from C64-lines")  # on/off flag
+parser.add_argument('-HC', action='store_true', help = "evaluates housekeeping data from HC-lines")  # on/off flag
+parser.add_argument('-H', action='store_true', help = "evaluates housekeeping data from H-lines")  # on/off flag
+
+parser.add_argument("-nameadd", type=str, nargs='+', help = "custom nameadd to add to filename")
+
+args = parser.parse_args()
+
+file = args.file
+if not Path(file).is_file():
+    print("no valid file: "+str(file))
+    sys.exit()
+
+filename = file.split(".")[0].split("/")[-1]
+
+if args.verbose == None:
+    verbosity = 1
+else:
+    verbosity = args.verbose
+
+
+def main():
+    
+    if args.Itime:
+        Do_Itime()
+        print('.Itime was created')
+               
+    if args.hist:
+        Do_hist()
+        print('.hist was created')
+
+    if args.HC:
+        Do_HC()
+        print('.hc was created')
+    
+    if args.H:
+        Do_H()
+        print('.h was created')
+        
+    if args.P:
+        Do_P()
+        print('.p was created')
+        
+    if args.C64:
+        Do_C64()
+        print('.c64 was created')
+            
+
+def Do_hist():
+    global mV; mV = 14000
+    global minV; minV = -100
+    global maxV; maxV = 3500
+    global resV; resV = 2*0.838214/2
+    global pol; pol = [1.88909231e-04, -6.84421799e-03, -3.28911747e+00,  3.24630271e+02]
+
+    global thr; global ch; global trig; trig=[]
+    ch,name,thr,u = channels()
+    chans = list(range(0, 18))
+    nameadd='hk'
+
+    energy = 'keV'
+    if args.mV: 
+        energy = 'mV'
+        nameadd = nameadd + '_mV'
+        u = [1.0]*18
+        
+    if args.trigger:
+        nameadd = nameadd + '_trigger'
+        for t in args.trigger:
+            trig.append(t) 
+            nameadd = nameadd + '_' + str(t)
+            
+    if args.xray:
+        nameadd = nameadd + '_xray'
+        
+    if not args.nameadd == None:
+        nameadd = args.nameadd[0]
+
+    global trigchan 
+    trigchan = []
+    for t in trig:
+        for i in range(len(name)):
+            if t == name[i]:
+                trigchan.append(i)
+
+    bin = int((maxV-minV)/resV)
+    
+    hist = np.zeros((bin+1,18+1))
+  
+    for i in range(0,bin+1):
+        hist[i][0]=minV + resV*i
+
+    with fileinput.input(files=(file), encoding="utf-8", errors='ignore') as f: 
+
+        for line in f:
+            l = line.split()
+            if l == []:
+                continue
+            if not 'EOF' in l and not 'X' in l and not 'xF' in l:
+                if l[0] == 'H':
+                    pass
+                    #T = adc_to_T(int(l[8]))
+                elif l[0] == 'EI' and len(l)>18:    
+                    if conditions(l,args):  
+                        for i in range(len(chans)):
+                            try: b = CmV(l,ch[chans[i]])*u[chans[i]] 
+                            except: b=maxV+10
+                            if b>=minV and b<=maxV:
+                                xb = int((b-minV)/resV)
+                                hist[xb,i+1] += 1                                                               
+                                
+    frame = pd.DataFrame(hist, columns = [energy] + name)
+    frame.to_csv('histograms/'+filename+nameadd+'.hist', sep = ' ', index = False)
+ 
+
+def Do_HC():
+    times, Vprims, Iprims, Ibiass, Vanas, Vdrvs, Idrvs, Tntcs, HVs = [], [], [], [], [], [], [], [], []
+# Vprim: Primary voltage
+# Iprim: Primary current provided
+# Ibias: Bias current (for the detectors)
+# Vana: 5 V for dac/adc (ana for analog?)
+# Vdrv: Voltage at input of  HV driver (Stephans resonator)
+# Idrv: Current at inpu  of HV driver (Stephans resonator)
+# Tntc: Temperature of NTC on Myrdins board
+# HV: High Voltage
+    R25 = 10e3
+    b = 3758
+    R1 = 33e3
+
+    with fileinput.input(files=(file), encoding="utf-8", errors = 'ignore') as f:
+          for line in f:
+              l = line.split()
+              if l == [] or 'EOF' in l or 'X' in l or  'xF' in l:
+                  continue
+              elif l[0] == 'HC':
+                  timestamp = int(l[1])
+                  if 1577836800 < timestamp < 2208985200:  # between 2020-01-01 and 2040-01-01
+                      time = datetime.fromtimestamp(float(l[1]))
+                      Vprim = (float(l[17]) / 0x10000) * 5.12 * 11.0
+                      Iprim = (float(l[5]) / 0x8000 - 1) * 0.64 * 2
+                      Ibias = (float(l[23]) / 0x8000 - 1) * 5.12 * 95.5
+                      Vana = (float(l[20]) / 0x10000) * 5.12 * 4.3
+                      Vdrv = (float(l[8]) / 0x10000) * 10.24
+                      Idrv = ((float(l[26]) / 0x10000) * 10.24 - Vdrv) / 100
+                      Vntc = (float(l[14]) / 0x10000) * 5.12
+                      Rntc = 33e3 * Vntc / (Vana - Vntc)
+                      Tntc = b / (np.log(Rntc / R25) + b / 298) - 273
+                      HV = (float(l[11]) / 0x8000 - 1) * 5.12 * 281
+
+                # Applying cuts for valid voltages and currents
+                  validV = all(0 <= v < 100 for v in [Vprim, Vana, Vdrv]) and (0 <= HV <= 1100) and (-20 < Tntc < 150)
+                  validA = all(0 <= c < 1000 for c in [Iprim, Ibias]) and -0.01 <= Idrv < 0.01
+
+                  if validV and validA:
+                      times.append(time)
+                      Vprims.append(Vprim)
+                      Iprims.append(Iprim)
+                      Ibiass.append(Ibias)
+                      Vanas.append(Vana)
+                      Vdrvs.append(Vdrv)
+                      Idrvs.append(Idrv)
+                      Tntcs.append(Tntc)
+                      HVs.append(HV)
+                      
+    frame = pd.DataFrame([times, Vprims, Iprims, Ibiass, Vanas, Vdrvs, Idrvs, Tntcs, HVs]).transpose()
+    frame.columns = ["time", "Vprim", "Iprim", "Ibias", "Vana", "Vdrv", "Idrv", "Tntc", "HV"]
+    frame.to_csv('hk_data/'+filename+'.hc', sep = ' ', index = False)
+ 
+#############################################################
+def degC(d, R1, R25, beta):
+    R = R1 * d / (1 - d)
+    T = beta / (np.log(R / R25) + beta / 298.0) - 273
+    return T
+
+def Do_H():
+    R25 = 10e3
+    b = 3758
+    R1 = 10e3
+    R1_8_12 = 33e3
+    T_ref = 65536
+    times, Ibiass, Vbiass, Iccs, Isss, Vccs, Vsss, Ts, Tavas, T8s, T12s = [], [], [], [], [], [], [], [], [],[],[]
+    ref = 3.3/65536
+    with fileinput.input(files=(file), encoding="utf-8", errors = 'ignore') as f:
+          for line in f:
+              l = line.split()
+              if l == [] or 'EOF' in l or 'X' in l or  'xF' in l:
+                  continue
+              elif l[0] == 'H':
+                  timestamp = int(l[1])
+                  if 1577836800 < timestamp < 2208985200:  # between 2020-01-01 and 2040-01-01
+                      time = datetime.fromtimestamp(float(l[1]))
+                      Ibias = float(l[4]) * ref * 97
+                      Vbias = -1 * float(l[10]) * ref * 50
+                      Icc = float(l[5]) * ref * 100
+                      Iss = -1 * float(l[6]) * ref * 100
+                      Vcc = float(l[12]) * ref * 2
+                      Vss = float(l[13]) * ref * 2
+
+                      T8_R1 = degC(float(l[7]) / 0x10000, R1_8_12, R25, b)  # Temp Sensor at position 8
+                      T12_R1 = degC(float(l[11]) / 0x10000, R1_8_12, R25, b)  # Temp Sensor at position 12
+
+                      T = float(l[8]) / T_ref  # Avas BGO
+                      T = R1 * T / (1 - T)
+                      T = b / (np.log(T/R25) + b/298.0) - 273
+                      Tava = adc_to_T(int(l[8]))
+
+                # Applying cuts for valid voltages and currents
+                  validV = all([-100 <= Vbias <= 0.0, 0 <= Vcc <= 100, 0 <= Vss <= 100, -20 <= T8_R1 <= 150, -20 <= T12_R1 <= 150, -20 <= T <= 150])
+                  validA = all([0 <= Ibias < 500, 0 <= Icc <= 50, -50 <= Iss <= 0])
+                  if validV and validA:
+                            times.append(time)
+                            Ibiass.append(Ibias)
+                            Vbiass.append(Vbias)
+                            Iccs.append(Icc)
+                            Isss.append(Iss)
+                            Vccs.append(Vcc)
+                            Vsss.append(Vss)
+                            Ts.append(T)
+                            Tavas.append(Tava)
+                            T8s.append(T8_R1)
+                            T12s.append(T12_R1)
+                            
+                            
+    frame = pd.DataFrame([times, Ibiass, Vbiass, Iccs, Isss, Vccs, Vsss, Ts, Tavas,T8s,T12s]).transpose()
+    frame.columns = ["time", "Ibias", "Vbias", "Icc", "Iss", "Vcc", "Vss", "T", "Tava","T8","T12"]
+    frame.to_csv('hk_data/'+filename+'.h', sep = ' ', index = False)
+    
+####################################################################
+
+def Do_P():
+    t1=[];t2=[];p1=[];p2=[];tstamps1=[];tstamps2=[]
+    time = None
+    
+    with fileinput.input(files=(file), encoding="utf-8", errors = 'ignore') as f:
+        for line in f:
+            l = line.split()
+            if not 'EOF' in l and not 'X' in l and not 'xF' in l and not l == []:
+                if l[0] == 'H':
+                    #T_BGO = adc_to_T(int(l[8]))
+                    timestamp = l[1]
+                    time = datetime.fromtimestamp(int(l[1]))
+                    #time.append(l[1])
+                    #temp.append(T_BGO)
+                    
+                elif l[0] == 'P' and time != None:
+                    hex_int = lambda x: int(x, 16)
+                    Tcal, pcal = senscalib([hex_int(l[3]),hex_int(l[4]),hex_int(l[5]),hex_int(l[6])], hex_int(l[8]), hex_int(l[7]))
+                    if l[1] == '0xba7e':
+                        process_sensor_data(time,Tcal, pcal, t1, p1, tstamps1)
+                    elif l[1] == '0xbafe':
+                        process_sensor_data(time,Tcal, pcal, t2, p2, tstamps2)
+                    
+    frame = pd.DataFrame([tstamps1,t1,p1,tstamps2,t2,p2]).transpose()
+    frame.columns = ["time1","T1","p1","time2","T2","p2"]
+    frame.to_csv('hk_data/'+filename+'.p', sep = ' ', index = False)
+    
+def Do_C64():
+    counts = [[] for i in range(18)]
+    times = []
+    time = None
+    ch,name,thr,u = channels()
+    
+    with fileinput.input(files=(file), encoding="utf-8", errors = 'ignore') as f:
+        for line in f:
+            l = line.split()
+            if not 'EOF' in l and not 'X' in l and not 'xF' in l and not l == []:
+                if l[0] == 'H':
+                    timestamp = l[1]
+                    time = datetime.fromtimestamp(int(l[1]))
+                    
+                elif l[0] == 'C64' and time != None:
+                    for i in range(18):
+                        counts[i].append(float(l[ch[i]+2])/10) #counts per 10s -> divide by 10 to get counts per second
+                    times.append(time) 
+                    
+    frame = pd.DataFrame([times]+counts).transpose()
+    frame.columns = ["time"]+name
+    frame.to_csv('hk_data/'+filename+'.c64', sep = ' ', index = False)
+ 
+
+ 
+def process_sensor_data(timestamp, Tcal, pcal, Ts, ps, stamps):
+    if -100 <= Tcal <= 100 and 0 <= pcal <= 1200:
+        Ts.append(Tcal)
+        ps.append(pcal)
+    elif -100 <= Tcal <= 100:
+        Ts.append(Tcal)
+        ps.append(None)
+    elif 0 <= pcal <= 1200:
+        Ts.append(None)
+        ps.append(pcal)
+    
+    stamps.append(timestamp)
+    
+def senscalib(word,T,p): #calibrates temperature and pressure values from MS5534C
+    # Convert calibration data into coefficients
+    C0 = int(word[0]) >> 1
+    C1 = ((word[2] & 0x3F) << 6) | (word[3] & 0x3F)
+    C2 = word[3] >> 6
+    C3 = word[2] >> 6
+    C4 = ((word[0] & 1) << 10) | (word[1] >> 6)
+    C5 = word[1] & 0x3F
+
+    # Calculate actual temperature
+    dT = T - (8 * C4 + 20224)
+    TEMP = 200 + dT * (C5 + 50) / 1024
+
+    # Calculate temperature compensated pressure
+    OFF = C1 * 4 + ((C3 - 512) * dT) / 4096
+    SENS = C0 + (C2 * dT) / 1024 + 24576
+    P = (SENS * (p - 7168)) / 16384 - OFF
+    P = P * 10 / 32 + 2500
+
+    # Second-order correction for non-linearity
+    if TEMP < 200:
+        T2 = 11 * (C5 + 24) * (200 - TEMP) ** 2 / 1048576
+        P2 = 3 * T2 * (P - 3500) / 16384
+    elif TEMP > 450:
+        T2 = 3 * (C5 + 24) * (450 - TEMP) ** 2 / 1048576
+        P2 = T2 * (P - 10000) / 8192
+    else:
+        T2 = P2 = 0
+    TEMP -= T2
+    P -= P2
+
+    return TEMP * 0.1, P * 0.1
+    
+################################################################
+
+def Do_Itime():
+    Itime = find_Itime()
+    fname = 'histograms/'+filename+ '.Itime'       
+    outfile = open(fname, 'w')
+    outfile.write(str(Itime))
+    outfile.close()
+
+def find_Itime():
+    Itime = 0
+    lasttime = 0
+    with fileinput.input(files=(file), encoding="utf-8", errors = 'ignore') 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 find_timentc():
+    time=[]
+    ntc=[]
+    with fileinput.input(files=(file), encoding="utf-8", errors = 'ignore') as f:
+        for line in f:
+            l = line.split() 
+            if not 'EOF' in l and l[0] == 'H':
+                T = l[8]
+                time.append(l[1])
+                ntc.append(T)
+    return time,ntc
+
+def find_timetemp():
+    time=[]
+    temp=[]
+    with fileinput.input(files=(file), encoding="utf-8", errors = 'ignore') as f:
+        for line in f:
+            l = line.split() 
+            if not 'EOF' in l and l[0] == 'H':
+                T = adc_to_T(int(l[8]))
+                time.append(l[1])
+                temp.append(T)
+    return time,temp
+
+########################################################################
+def channels():
+    ch = np.arange(18)
+    name=[None]*18
+    thr=[12]*18
+    u = [1.0]*18
+    s = [1.0]*18
+   
+    name[0] = "A1H";  ch[0] = 13;  thr[0] = 10;  s[0] = 15.0;  u[0] = 4.55
+    name[1] = "A1L";  ch[1] = 12;  thr[1] = 10;  s[1] = 1.0;   u[1] = 1.0
+    name[2] = "A2H";  ch[2] = 10;  thr[2] = 15;  s[2] = 15.0;  u[2] = 4.55
+    name[3] = "A2L";  ch[3] = 9;   thr[3] = 10;  s[3] = 1.0;   u[3] = 1.0                    
+    name[4] = "B";    ch[4] = 17;  thr[4] = 50;  s[4] = 15.0;  u[4] = 1.0                    
+    name[5] = "C1H";  ch[5] = 7;   thr[5] = 10;  s[5] = 15.0;  u[5] = 4.55                    
+    name[6] = "C1L";  ch[6] = 6;   thr[6] = 10;  s[6] = 1.0;   u[6] = 1.0                   
+    name[7] = "C2H";  ch[7] = 4;   thr[7] = 15;  s[7] = 15.0;  u[7] = 4.55                  
+    name[8] = "C2L";  ch[8] = 3;   thr[8] = 10;  s[8] = 1.0;   u[8] = 1.0                   
+    name[9] = "D1H";  ch[9] = 16;  thr[9] = 20;  s[9] = 15.0;  u[9] = 1.577                   
+    name[10] = "D1L"; ch[10] = 15; thr[10] = 20; s[10] = 1.0;  u[10] = 1.0
+    name[11] = "D2H"; ch[11] = 14; thr[11] = 20; s[11] = 15.0; u[11] = 1.577
+    name[12] = "D2L"; ch[12] = 11; thr[12] = 20; s[12] = 1.0;  u[12] = 1.0
+    name[13] = "EH";  ch[13] = 1;  thr[13] = 10; s[13] = 15.0; u[13] = 4.55
+    name[14] = "EL";  ch[14] = 0;  thr[14] = 10; s[14] = 1.0;  u[14] = 1.0
+    name[15] = "E1";  ch[15] = 8;  thr[15] = 8;  s[15] = 2.2;  u[15] = 11.32
+    name[16] = "E2";  ch[16] = 5;  thr[16] = 8;  s[16] = 2.2;  u[16] = 11.32
+    name[17] = "E3";  ch[17] = 2;  thr[17] = 8;  s[17] = 2.2;  u[17] = 11.32
+   
+    return ch, name, thr, u 
+
+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 #bis 2024-05-28
+    #R25 = 9820
+    #R1 = 30700
+    b = 3940
+    #dac = 62690
+    dac = 62697 #new value from Stephan (see matrix, AHBGO 28.05.24) #bis 2024-05-28
+    #dac = 58380
+    T = 0
+    #----------------------Daten für CHAOS von Hannes
+    R25 = 10e3
+    b = 3758
+    R1 = 33e3
+
+    alpha = np.log(R1/R25 *(dac/(dac-adc)-1))
+    T = 1/(alpha/b+1/298.15) - 273.15
+    return T 
+
+###############################################################
+
+def conditions(l,args):
+    try: 
+        if not args.trigger or isTrigger(l):
+            if not args.xray or not isBGO(l):
+                return True
+        else: return False
+    except: return False
+    
+def isTrigger(l):
+    for t in trigchan:
+        if CmV(l,ch[t]) <= thr[t]: return False
+    return True
+
+def isBGO(l):
+    if CmV(l,ch[9]) > thr[9] and CmV(l,ch[11]) > thr[11]: return True
+    return False
+
+main()

CHAOSx.py

@@ -13,15 +13,12 @@ parser.add_argument("file", type = str, help="path/file were data is stored in")
 parser.add_argument('-Itime', action='store_true', help = "if chosen, Itime is calculatet and printet to filename.Itime")  # on/off flag
 
 parser.add_argument('-hist', action='store_true', help = "makes histogramm of all channels using EI-lines")  # on/off flag
-parser.add_argument("-mV", action='store_true', help = "needs option -hist; if chosen, data will be kept in mV changed to keV otherwise")  # on/off flag
+parser.add_argument('-tdhist', type=str, nargs='+', help = "need two arguments: channelx, channely -> makes 2D histogramm for channelx and channely")  # on/off flag
+
+parser.add_argument("-mV", action='store_true', help = "needs option -hist or -tdhist; if chosen, data will be kept in mV changed to keV otherwise")  # on/off flag
 parser.add_argument("-trigger", type=str, nargs='+', help = "needs option -EI; letters of Triggerchannels to chose")
 parser.add_argument("-xray", action='store_true', help = "needs option -EI; only select xray data")  # on/off flag
 
-parser.add_argument('-P', action='store_true', help = "evaluates pressure and temperature data from P-lines")  # on/off flag
-parser.add_argument('-C64', action='store_true', help = "evaluates counts per second from C64-lines")  # on/off flag
-parser.add_argument('-HC', action='store_true', help = "evaluates housekeeping data from HC-lines")  # on/off flag
-parser.add_argument('-H', action='store_true', help = "evaluates housekeeping data from H-lines")  # on/off flag
-
 parser.add_argument("-nameadd", type=str, nargs='+', help = "custom nameadd to add to filename")
 
 args = parser.parse_args()
@@ -49,21 +46,9 @@ def main():
         Do_hist()
         print('.hist was created')
 
-    if args.HC:
-        Do_HC()
-        print('.hc was created')
-    
-    if args.H:
-        Do_H()
-        print('.h was created')
-        
-    if args.P:
-        Do_P()
-        print('.p was created')
-        
-    if args.C64:
-        Do_C64()
-        print('.c64 was created')
+    if args.tdhist:
+        Do_2dhist(args.tdhist)
+        print('.2dhist was created')
             
 
 def Do_hist():
@@ -114,7 +99,7 @@ def Do_hist():
 
         for line in f:
             l = line.split()
-            if l == []:
+            if l == [] or len(l)<18:
                 continue
             if not 'EOF' in l and not 'X' in l and not 'xF' in l:
                 if l[0] == 'H':
@@ -131,220 +116,85 @@ def Do_hist():
     frame = pd.DataFrame(hist, columns = [energy] + name)
     frame.to_csv('histograms/'+filename+nameadd+'.hist', sep = ' ', index = False)
  
+def Do_2dhist(chans):
+    global mV; mV = 14000
+    global minV; minV = -10
+    global maxV; maxV = 2500
+    global resV; resV = 5*0.838214
+    global pol; pol = [1.88909231e-04, -6.84421799e-03, -3.28911747e+00,  3.24630271e+02]
 
-def Do_HC():
-    times, Vprims, Iprims, Ibiass, Vanas, Vdrvs, Idrvs, Tntcs, HVs = [], [], [], [], [], [], [], [], []
-# Vprim: Primary voltage
-# Iprim: Primary current provided
-# Ibias: Bias current (for the detectors)
-# Vana: 5 V for dac/adc (ana for analog?)
-# Vdrv: Voltage at input of  HV driver (Stephans resonator)
-# Idrv: Current at inpu  of HV driver (Stephans resonator)
-# Tntc: Temperature of NTC on Myrdins board
-# HV: High Voltage
-    R25 = 10e3
-    b = 3758
-    R1 = 33e3
-
-    with fileinput.input(files=(file), encoding="utf-8", errors = 'ignore') as f:
-          for line in f:
-              l = line.split()
-              if l == [] or 'EOF' in l or 'X' in l or  'xF' in l:
-                  continue
-              elif l[0] == 'HC':
-                  timestamp = int(l[1])
-                  if 1577836800 < timestamp < 2208985200:  # between 2020-01-01 and 2040-01-01
-                      time = datetime.fromtimestamp(float(l[1]))
-                      Vprim = (float(l[17]) / 0x10000) * 5.12 * 11.0
-                      Iprim = (float(l[5]) / 0x8000 - 1) * 0.64 * 2
-                      Ibias = (float(l[23]) / 0x8000 - 1) * 5.12 * 95.5
-                      Vana = (float(l[20]) / 0x10000) * 5.12 * 4.3
-                      Vdrv = (float(l[8]) / 0x10000) * 10.24
-                      Idrv = ((float(l[26]) / 0x10000) * 10.24 - Vdrv) / 100
-                      Vntc = (float(l[14]) / 0x10000) * 5.12
-                      Rntc = 33e3 * Vntc / (Vana - Vntc)
-                      Tntc = b / (np.log(Rntc / R25) + b / 298) - 273
-                      HV = (float(l[11]) / 0x8000 - 1) * 5.12 * 281
-
-                # Applying cuts for valid voltages and currents
-                  validV = all(0 <= v < 100 for v in [Vprim, Vana, Vdrv]) and (0 <= HV <= 1100) and (-20 < Tntc < 150)
-                  validA = all(0 <= c < 1000 for c in [Iprim, Ibias]) and -0.01 <= Idrv < 0.01
-
-                  if validV and validA:
-                      times.append(time)
-                      Vprims.append(Vprim)
-                      Iprims.append(Iprim)
-                      Ibiass.append(Ibias)
-                      Vanas.append(Vana)
-                      Vdrvs.append(Vdrv)
-                      Idrvs.append(Idrv)
-                      Tntcs.append(Tntc)
-                      HVs.append(HV)
-                      
-    frame = pd.DataFrame([times, Vprims, Iprims, Ibiass, Vanas, Vdrvs, Idrvs, Tntcs, HVs]).transpose()
-    frame.columns = ["time", "Vprim", "Iprim", "Ibias", "Vana", "Vdrv", "Idrv", "Tntc", "HV"]
-    frame.to_csv('hk_data/'+filename+'.hc', sep = ' ', index = False)
- 
-#############################################################
-def degC(d, R1, R25, beta):
-    R = R1 * d / (1 - d)
-    T = beta / (np.log(R / R25) + beta / 298.0) - 273
-    return T
-
-def Do_H():
-    R25 = 10e3
-    b = 3758
-    R1 = 10e3
-    R1_8_12 = 33e3
-    T_ref = 65536
-    times, Ibiass, Vbiass, Iccs, Isss, Vccs, Vsss, Ts, Tavas, T8s, T12s = [], [], [], [], [], [], [], [], [],[],[]
-    ref = 3.3/65536
-    with fileinput.input(files=(file), encoding="utf-8", errors = 'ignore') as f:
-          for line in f:
-              l = line.split()
-              if l == [] or 'EOF' in l or 'X' in l or  'xF' in l:
-                  continue
-              elif l[0] == 'H':
-                  timestamp = int(l[1])
-                  if 1577836800 < timestamp < 2208985200:  # between 2020-01-01 and 2040-01-01
-                      time = datetime.fromtimestamp(float(l[1]))
-                      Ibias = float(l[4]) * ref * 97
-                      Vbias = -1 * float(l[10]) * ref * 50
-                      Icc = float(l[5]) * ref * 100
-                      Iss = -1 * float(l[6]) * ref * 100
-                      Vcc = float(l[12]) * ref * 2
-                      Vss = float(l[13]) * ref * 2
-
-                      T8_R1 = degC(float(l[7]) / 0x10000, R1_8_12, R25, b)  # Temp Sensor at position 8
-                      T12_R1 = degC(float(l[11]) / 0x10000, R1_8_12, R25, b)  # Temp Sensor at position 12
-
-                      T = float(l[8]) / T_ref  # Avas BGO
-                      T = R1 * T / (1 - T)
-                      T = b / (np.log(T/R25) + b/298.0) - 273
-                      Tava = adc_to_T(int(l[8]))
-
-                # Applying cuts for valid voltages and currents
-                  validV = all([-100 <= Vbias <= 0.0, 0 <= Vcc <= 100, 0 <= Vss <= 100, -20 <= T8_R1 <= 150, -20 <= T12_R1 <= 150, -20 <= T <= 150])
-                  validA = all([0 <= Ibias < 500, 0 <= Icc <= 50, -50 <= Iss <= 0])
-                  if validV and validA:
-                            times.append(time)
-                            Ibiass.append(Ibias)
-                            Vbiass.append(Vbias)
-                            Iccs.append(Icc)
-                            Isss.append(Iss)
-                            Vccs.append(Vcc)
-                            Vsss.append(Vss)
-                            Ts.append(T)
-                            Tavas.append(Tava)
-                            T8s.append(T8_R1)
-                            T12s.append(T12_R1)
-                            
-                            
-    frame = pd.DataFrame([times, Ibiass, Vbiass, Iccs, Isss, Vccs, Vsss, Ts, Tavas,T8s,T12s]).transpose()
-    frame.columns = ["time", "Ibias", "Vbias", "Icc", "Iss", "Vcc", "Vss", "T", "Tava","T8","T12"]
-    frame.to_csv('hk_data/'+filename+'.h', sep = ' ', index = False)
-    
-####################################################################
-
-def Do_P():
-    t1=[];t2=[];p1=[];p2=[];tstamps1=[];tstamps2=[]
-    time = None
-    
-    with fileinput.input(files=(file), encoding="utf-8", errors = 'ignore') as f:
-        for line in f:
-            l = line.split()
-            if not 'EOF' in l and not 'X' in l and not 'xF' in l and not l == []:
-                if l[0] == 'H':
-                    #T_BGO = adc_to_T(int(l[8]))
-                    timestamp = l[1]
-                    time = datetime.fromtimestamp(int(l[1]))
-                    #time.append(l[1])
-                    #temp.append(T_BGO)
-                    
-                elif l[0] == 'P' and time != None:
-                    hex_int = lambda x: int(x, 16)
-                    Tcal, pcal = senscalib([hex_int(l[3]),hex_int(l[4]),hex_int(l[5]),hex_int(l[6])], hex_int(l[8]), hex_int(l[7]))
-                    if l[1] == '0xba7e':
-                        process_sensor_data(time,Tcal, pcal, t1, p1, tstamps1)
-                    elif l[1] == '0xbafe':
-                        process_sensor_data(time,Tcal, pcal, t2, p2, tstamps2)
-                    
-    frame = pd.DataFrame([tstamps1,t1,p1,tstamps2,t2,p2]).transpose()
-    frame.columns = ["time1","T1","p1","time2","T2","p2"]
-    frame.to_csv('hk_data/'+filename+'.p', sep = ' ', index = False)
-    
-def Do_C64():
-    counts = [[] for i in range(18)]
-    times = []
-    time = None
+    global thr; global ch; global trig; trig=[]
     ch,name,thr,u = channels()
+    x = name.index(chans[0])
+    y = name.index(chans[1])
+
+    nameadd='x'
+
+    energy = 'keV'
+    if args.mV: 
+        energy = 'mV'
+        nameadd = nameadd + '_mV'
+        u = [1.0]*18
+        
+    if args.trigger:
+        nameadd = nameadd + '_trigger'
+        for t in args.trigger:
+            trig.append(t) 
+            nameadd = nameadd + '_' + str(t)
+            
+    if args.xray:
+        nameadd = nameadd + '_xray'
+        
+    if not args.nameadd == None:
+        nameadd = args.nameadd[0]
+
+    global trigchan 
+    trigchan = []
+    for t in trig:
+        for i in range(len(name)):
+            if t == name[i]:
+                trigchan.append(i)
+
+    bin = int((maxV-minV)/resV)
     
-    with fileinput.input(files=(file), encoding="utf-8", errors = 'ignore') as f:
+    #hist = np.zeros(((bin+1)*(bin+1),3))
+    hist = np.zeros((bin+1,bin+2))
+    edges = []
+    for i in range(0,bin+1):
+        # for j in range(bin+1):
+        #     hist[i*(bin+1)+j][0] = minV + resV*i
+        #     hist[i*(bin+1)+j][1] = minV + resV*j
+        edge = minV + resV*i
+        hist[i][0] = edge
+        edges.append(edge)
+            
+
+    with fileinput.input(files=(file), encoding="utf-8", errors='ignore') as f: 
+
         for line in f:
             l = line.split()
-            if not 'EOF' in l and not 'X' in l and not 'xF' in l and not l == []:
+            if l == []:
+                continue
+            if not 'EOF' in l and not 'X' in l and not 'xF' in l:
                 if l[0] == 'H':
-                    timestamp = l[1]
-                    time = datetime.fromtimestamp(int(l[1]))
-                    
-                elif l[0] == 'C64' and time != None:
-                    for i in range(18):
-                        counts[i].append(float(l[ch[i]+2])/10) #counts per 10s -> divide by 10 to get counts per second
-                    times.append(time) 
-                    
-    frame = pd.DataFrame([times]+counts).transpose()
-    frame.columns = ["time"]+name
-    frame.to_csv('hk_data/'+filename+'.c64', sep = ' ', index = False)
- 
+                    pass
+                    #T = adc_to_T(int(l[8]))
+                elif l[0] == 'EI' and len(l)>3*18:    
+                    if conditions(l,args):  
+                        bx = CmV(l,ch[x])*u[x] 
+                        by = CmV(l,ch[y])*u[y] 
 
- 
-def process_sensor_data(timestamp, Tcal, pcal, Ts, ps, stamps):
-    if -100 <= Tcal <= 100 and 0 <= pcal <= 1200:
-        Ts.append(Tcal)
-        ps.append(pcal)
-    elif -100 <= Tcal <= 100:
-        Ts.append(Tcal)
-        ps.append(None)
-    elif 0 <= pcal <= 1200:
-        Ts.append(None)
-        ps.append(pcal)
-    
-    stamps.append(timestamp)
-    
-def senscalib(word,T,p): #calibrates temperature and pressure values from MS5534C
-    # Convert calibration data into coefficients
-    C0 = int(word[0]) >> 1
-    C1 = ((word[2] & 0x3F) << 6) | (word[3] & 0x3F)
-    C2 = word[3] >> 6
-    C3 = word[2] >> 6
-    C4 = ((word[0] & 1) << 10) | (word[1] >> 6)
-    C5 = word[1] & 0x3F
-
-    # Calculate actual temperature
-    dT = T - (8 * C4 + 20224)
-    TEMP = 200 + dT * (C5 + 50) / 1024
-
-    # Calculate temperature compensated pressure
-    OFF = C1 * 4 + ((C3 - 512) * dT) / 4096
-    SENS = C0 + (C2 * dT) / 1024 + 24576
-    P = (SENS * (p - 7168)) / 16384 - OFF
-    P = P * 10 / 32 + 2500
-
-    # Second-order correction for non-linearity
-    if TEMP < 200:
-        T2 = 11 * (C5 + 24) * (200 - TEMP) ** 2 / 1048576
-        P2 = 3 * T2 * (P - 3500) / 16384
-    elif TEMP > 450:
-        T2 = 3 * (C5 + 24) * (450 - TEMP) ** 2 / 1048576
-        P2 = T2 * (P - 10000) / 8192
-    else:
-        T2 = P2 = 0
-    TEMP -= T2
-    P -= P2
-
-    return TEMP * 0.1, P * 0.1
-    
-################################################################
+                        if bx>=minV and bx<=maxV and by>=minV and by<=maxV:
+                            # bxx = int((bx-minV)/resV)
+                            # byy = int((by-minV)/resV)
+                            # hist[bxx*(bin+1)+byy,2] += 1     
+                            bxx = int((bx-minV)/resV)
+                            byy = int((by-minV)/resV)+1
+                            hist[bxx,byy] += 1                                                        
+                                
+    #frame = pd.DataFrame(hist, columns = [name[x],name[y],'counts'])
+    frame = pd.DataFrame(hist, columns = ['keV']+edges)
+    frame.to_csv('histograms/'+filename+nameadd+'_'+name[x]+'vs'+name[y]+'.2dhist', sep = ' ', index = False)
 
 def Do_Itime():
     Itime = find_Itime()
@@ -448,14 +298,12 @@ def adc_to_T(adc):
 
 ###############################################################
 
-def conditions(l,args):
+def conditions(l,args): 
     if not args.trigger or isTrigger(l):
-        try:
-            if not args.xray or not isBGO(l):
-                return True
-        except: return False
+        if not args.xray or not isBGO(l):
+            return True
     else: return False
-    
+
 def isTrigger(l):
     for t in trigchan:
         if CmV(l,ch[t]) <= thr[t]: return False

fit_edges.py

@@ -25,7 +25,7 @@ def main():
     file2 = "histograms/2024-07-22-Instrument-22V-1x_mV_xray.hist"
     #mV A1H A1L A2H A2L B C1H C1L C2H C2L D1H D1L D2H D2L EH EL E1 E2 E3
     
-    #data1, Itime1 = read_data(file1)
+    data1, Itime1 = read_data(file1)
     data2, Itime2 = read_data(file2)
     
     channels1 = ["A1H","A2H","C1H","C2H","EH"]
@@ -56,14 +56,25 @@ def main():
     
     plt.figure(figsize=(17,12))
     #fit_and_plot_data(data1,Itime1,channels1,"",a1,b1,startparams1)
-    fit_and_plot_data(data2,Itime2,channels2,"",a1,b1,startparams1)
+    #fit_and_plot_data(data2,Itime2,channels2,"",a1,b1,startparams1)
     
-    tit = "Fit of Compton edges for E123 (bi-source below; 07-22-Instrument)"
+###############################################################################
+
+    file = "histograms/2024-08-27-testing-2x_trigger_D1H_D2H.hist"
+    #mV A1H A1L A2H A2L B C1H C1L C2H C2L D1H D1L D2H D2L EH EL E1 E2 E3
+    
+    data, Itime = read_data(file)
+
+    channels = ['D1','D2']
+    
+    lowgain(data, Itime, channels)
+    
+    tit = "High Gain vs low gain"
     layout(tit)
     
-    #plt.show()
-    plt.savefig('plots/comptonE123.pdf')
-    plt.savefig('plots/comptonE123.png')
+    plt.show()
+    #plt.savefig('plots/comptonE123.pdf')
+    #plt.savefig('plots/comptonE123.png')
 
 
 def read_data(file):
@@ -102,6 +113,12 @@ def fit_and_plot_data(data, Itime, channels,add, a, b, startparams):
             plt.step(data['mV'], data[ch], label=ch)
             #plt.step(data['mV'], data[ch], label=ch+' ; $u_C$ = ('+ str(round(par[2]/c,4)) + '$\pm$' +str(round(perr[2]/c+0.00005,4))+')mV/keV')
             #plt.plot(data['mV'], y(data['mV']), 'k--') 
+            
+def lowgain(data, Itime, channels):
+    resV = 0.838214/2
+    for ch in channels:
+        plt.plot(np.multiply(data[ch+'H'], 3600/resV/Itime),np.multiply(data[ch+'L'], 3600/resV/Itime), label = ch)
+    
     
 def layout(tit): 
     fs = 15
@@ -111,9 +128,14 @@ def layout(tit):
     ylim = [0.3,30]
     xlim = [10,60] #E123
     ylim = [50,6000]
+    
+    #############
+    xlim=[0,300]
+    ylim=xlim
   
-    ylab = 'counts [1/h/mV]'
+    ylab = 'Low Gain counts [1/h/mV]'
     xlab = 'signal [mV]'
+    xlab='High Gain counts [1/h/mV]'
 
     plt.legend(fontsize = fs)
     plt.xticks(fontsize=fs)

histplot2d.py

@@ -0,0 +1,144 @@
+'''
+Autor: Ava Pohley
+'''
+
+import argparse
+import sys
+import matplotlib.pyplot as plt
+import pandas as pd
+import numpy as np
+import os.path
+from matplotlib import cm
+
+
+
+def main():
+    
+    parser = argparse.ArgumentParser(description="This script plots histogramms of given data ending on '.2dhist'")
+    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("filepath", type = str, help="path/file were data is stored in")
+    parser.add_argument("-t", "--title", type = str, help="title to use in plot; if not given a standard title will be chosen")
+    parser.add_argument("-e", "--export", type = str, nargs='+', help = "first: filepath were figure shall be exported to (if f, filepath of histdata will be chosen with filepath_plot) \n second:format (png,pdf,...)")
+    parser.add_argument('-r', '--raw', action='store_true', help = "just prints data without considering itime")  # on/off flag
+    parser.add_argument("-l", "--limits", type=str, nargs='+',help = "limits of plot in order xmin,xmax, ymin,ymax,")
+
+    args = parser.parse_args()
+
+    if args.verbose == None:
+        verbosity = 1
+    else:
+        verbosity = args.verbose
+    
+    if verbosity == 2:
+        sys.stderr.write("--verbostiy = '" + str(args.verbose)+"'\n")
+        sys.stderr.write("--filepath = '" + str(args.filepath)+"'\n")
+        sys.stderr.write("--title = '" + str(args.title)+"'\n")
+        sys.stderr.write("--export = '" + str(args.export) +"'\n")
+        sys.stderr.write("--raw = '" + str(args.raw) +"'\n")
+        sys.stderr.write("--limits = '" + str(args.limits) +"'\n")
+
+   
+    file = args.filepath.split('.')[0]
+    
+    end = '.2dhist'
+    energy = 'keV'
+    raw = args.raw
+    
+    if args.title == None:
+        tit = 'Histplot2D, '+"'"+file.split('/')[-1]+"' "
+    else:
+        tit = args.title       
+
+
+    x,y,c, data, Itime = read_data(file, end, raw)
+   
+    if Itime == None:
+        raw = True
+    
+    plt.figure(figsize=(15,10))
+
+    cbar = print_data(x,y,c,data,Itime, energy)
+
+    layout(x,y,cbar,tit,raw,args.limits,energy)
+
+
+    if args.export == None:
+        plt.show()
+    else:
+        if args.export[0] == 'f':
+            name = 'plots/plot_' + file.split('/')[-1]
+        else:
+            name = args.export[0]
+        if len(args.export)>1:
+                end = '.' + args.export[1]
+        else:
+                end=''
+        plt.savefig(name+end)
+
+
+def read_data(file, end, raw):
+
+    x,y,c = pd.read_csv(file+end, sep = ' ', nrows=0)
+    data = pd.read_csv(file+end, sep='\s+',skiprows=0)
+    
+    f = file + '.Itime'
+    Itime = None
+    if not raw == True:
+        while len(f) > 0:
+            if os.path.isfile(f+'.Itime'):
+                Itime = int(pd.read_csv(f+'.Itime', sep='\s+', skiprows=0, names = ['time'])['time'][0])
+                f=''
+            else:
+                f=f[:-1]
+
+    return x,y,c, data, Itime
+
+def print_data(x,y,c,data, Itime,energy): 
+    resV = data[y][1]-data[y][0]
+    
+    cmap = cm.get_cmap('jet', 256)
+    cmap.set_under('white')
+
+    im = plt.pcolormesh([data[x], data[y]], data[c], shading='auto', vmin=0.6)
+
+    #scatter = plt.scatter(data[x], data[y], c=data[c], vmin = 0.1) 
+    cbar = plt.colorbar(im, orientation = 'vertical', shrink = 0.95, pad = 0.1)
+   
+    return cbar
+        
+              
+
+def layout(x,y,cbar,tit,raw, limits,energy):
+
+    # energy = 'mV'
+    # if keV:
+    #     energy = 'keV'
+
+    fs = 15
+    lim = [None, None]
+
+    if not limits == None:
+        lim[0] = float(limits[0])
+        lim[1] = float(limits[1])
+    
+    if lim[0] == None:
+        lim[0] = -10
+    if lim[1] == None:
+        lim[1] = 350
+
+    ylab = y + 'in ' + energy+'$_{Si}$'
+    xlab = x + 'in ' + energy+'$_{Si}$'
+    
+    plt.xticks(fontsize=fs)
+    plt.yticks(fontsize=fs)
+    plt.title(tit, fontsize = fs+10, pad = 15)
+
+    plt.ylabel(ylab, fontsize = fs+5)
+    plt.xlabel(xlab, fontsize=fs+5)
+    plt.ylim(lim[0],lim[1])
+    plt.xlim(lim[0],lim[1])
+    
+    cbar.set_label('counts', fontsize = fs+5)
+    cbar.ax.tick_params(labelsize=fs)
+
+main()    

test2.py

@@ -0,0 +1,59 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Sep 23 10:40:21 2024
+
+@author: ava
+"""
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from matplotlib import cm
+
+
+minV = -2
+maxV = 4
+resV = 2
+bin = int((maxV-minV)/resV)
+
+hist = np.zeros(((bin+1)*(bin+1),3))
+columns = ['x','y','counts']
+for i in range(0,bin+1):
+    for j in range(bin+1):
+        hist[i*(bin+1)+j][0] = minV + resV*i
+        hist[i*(bin+1)+j][1] = minV + resV*j
+        
+for pair in [[2,4],[-2,2],[0,0],[0,0],[2,4],[0,0]]:    
+    bx,by=pair
+    bxx = int((bx-minV)/resV)
+    byy = int((by-minV)/resV)
+
+    hist[bxx*(bin+1)+byy,2] += 1  
+    
+print(hist)
+
+frame = pd.DataFrame(hist, columns = ['x_keV','y_keV','counts'])
+frame.to_csv('histograms/test2.2dhist', sep = ' ', index = False)
+
+# edges=[]
+#x,y,c = pd.read_csv('histograms/2024-08-27-testing-2X_D1HvsD2H.2dhist', sep = ' ', nrows=0)
+
+data = pd.read_csv('histograms/2024-08-27-testing-2X_D1HvsD2H.2dhist', sep='\s+',skiprows=0)
+
+c = np.array(data)[:,1:]
+print(c)
+# print(x)
+# print(y)
+# print(c)
+
+fig, ax = plt.subplots()
+cmap = cm.get_cmap()
+cmap.set_under('white')
+
+cb = plt.pcolormesh(data['keV'],data['keV'],c, cmap = cmap, shading='auto', vmin=0.1)
+#scatter = ax.scatter(data[x], data[y], c=data[c])
+# #plt.imshow(f[:,1:],interpolation='nearest',origin='upper')
+plt.colorbar(cb)
+# plt.hist2d(data['keV'],data['keV'],weights=counts)
+
+plt.show()