Entry 32106

# This script is for subtracting natural thorium from spectrum
# Calculation of half lives is also included

Directory = "C:\\HYPC\\Spectra\\"

# List of files

file_sam = ["a2th1p2","a2th2p2","a2th3p2",\
            "a2th4p2","a2th5p2","a2th6p2",\
            "a2th7p2","a2th8p2","a2th9p2"]
file_nat = "athnat.fon"

for k in range(0,len(file_sam)):
    
    # Open files
    filename_sam = Directory + file_sam[k] + ".pkl"
    filename_nat = Directory + file_nat
    sam = open(filename_sam, "r")
    nat = open(filename_nat, "r")
    # Read file in line by line
    sample = sam.readlines()
    natural = nat.readlines()
    # Close output file
    sam.close()
    nat.close()

    # Define variables for peak properties
    E_sam = []; E_sam_unc = []; FWHM_sam = []; i_sam = []; i_sam_rel = []
    E_nat = []; E_nat_unc = []; i_nat = []; i_nat_abs = []


    # Read data from sample
    for s in sample[14:len(sample)]:
        E_sam.append(float(s[28:37]))
        E_sam_unc.append(float(s[40:46]))
        FWHM_sam.append(float(s[48:53]))
        i_sam.append(float(s[61:72]))
        i_sam_rel.append(float(s[75:80]))

    # Read data from Natural Thorium
    for s in natural[4:len(natural)]:
        E_nat.append(float(s[1:10]))
        E_nat_unc.append(float(s[12:18]))
        i_nat.append(float(s[23:31]))
        i_nat_abs.append(float(s[36:41]))

    # Match the peak energies

    match = []
    nomatch = []
    areas_1 = []; areas_2 = []; areas_3 = []; areas_4 = [] #sam, sam_unc, nat, nat_unc
    vars()[file_sam[k]] = []
    #vars()[file_sam[k]+"_en"] = []
    #vars()[file_sam[k]+"_en_unc"] = []
    #vars()[file_sam[k]+"_area"] = []
    #vars()[file_sam[k]+"_area_unc"] = []

    j = 0
    for s in E_sam:
        for i in range(0,len(E_nat)):
            if (E_nat[i]-E_nat_unc[i]-E_sam_unc[j]-(0.5*FWHM_sam[j])) \
                     < s < \
               (E_nat[i]+E_nat_unc[i]+E_sam_unc[j]+(0.5*FWHM_sam[j])):
                match.append(s)
                #print "match  ", s
                areas_1.append(i_sam[j]); areas_2.append(i_sam_rel[j])
                areas_3.append(i_nat[i]); areas_4.append(i_nat_abs[i])
        if s not in match:
            nomatch.append(s)
            vars()[file_sam[k]].append([E_sam[j],E_sam_unc[j],\
                                        i_sam[j],i_sam_rel[j],k+1])
            #vars()[file_sam[k]+"_en"].append(E_sam[j])
            #vars()[file_sam[k]+"_en_unc"].append(E_sam_unc[j])
            #vars()[file_sam[k]+"_area"].append(i_sam[j])
            #vars()[file_sam[k]+"_area_unc"].append(i_sam_rel[j])
            #print "----nomatch  ", s
        j = j + 1
    #print len(vars()[file_sam[k]])
    #print (vars()[file_sam[k]])
    #print (sample[6][11:-5])
    vars()["livetime_"+file_sam[k]] = float(sample[6][11:-5])

livetime_nat = 143620.6

# Combine lists together

combined = []

for k in range(0,len(file_sam)):
    for s in vars()[file_sam[k]]:
        combined.append(s)
#print combined

# Sort all lists

sortedc = sorted(combined)

# Create new lists based on peak energy

i = 0 # position in sortedc
j = 1 # peak number

# If peak within boundaries of before peak then append to list
# If peak disagrees with previous peak then create new list
# Fifth index will be added and will be spectrum number

for i in range(0,len(sortedc)):
    if i == 0:
        vars()["Peak_"+str(j)] = []
        vars()["Peak_"+str(j)].append(sortedc[i])
    elif (sortedc[i][0]-sortedc[i-1][1]-sortedc[i][1])<\
         sortedc[i-1][0]<\
         sortedc[i][0]+sortedc[i-1][1]+sortedc[i][1]:
        vars()["Peak_"+str(j)].append(sortedc[i])
    else:
        if len(vars()["Peak_"+str(j)])<2:
            del vars()["Peak_"+str(j)]
            j = j-1
        
        #print j
        j = j+1
        vars()["Peak_"+str(j)] = []

for k in range(1,j):

    for s in (vars()["Peak_"+str(k)]):
        if s[4] == 1:
            s[2] = s[2]/livetime_a2th1p2
            s[3] = s[2]*s[3]/100
        elif s[4] == 2:
            s[2] = s[2]/livetime_a2th2p2
            s[3] = s[2]*s[3]/100
        elif s[4] == 3:
            s[2] = s[2]/livetime_a2th3p2
            s[3] = s[2]*s[3]/100
        elif s[4] == 4:
            s[2] = s[2]/livetime_a2th4p2
            s[3] = s[2]*s[3]/100
        elif s[4] == 5:
            s[2] = s[2]/livetime_a2th5p2
            s[3] = s[2]*s[3]/100
        elif s[4] == 6:
            s[2] = s[2]/livetime_a2th6p2
            s[3] = s[2]*s[3]/100
        elif s[4] == 7:
            s[2] = s[2]/livetime_a2th7p2
            s[3] = s[2]*s[3]/100
        elif s[4] == 8:
            s[2] = s[2]/livetime_a2th8p2
            s[3] = s[2]*s[3]/100
        elif s[4] == 9:
            s[2] = s[2]/livetime_a2th9p2
            s[3] = s[2]*s[3]/100

    #print k
    vars()["Peak_"+str(k)] = sorted(vars()["Peak_"+str(k)], key=lambda s: s[4])
    #print (vars()["Peak_"+str(k)])

# Plot decay curves

x = []
y = []
yerr = []

for s in Peak_1:
    y.append(s[2])
    yerr.append(s[3])
    if s[4] == 1:
        x.append( 4741 )
    if s[4] == 2:
        x.append( 10386 )
    if s[4] == 3:
        x.append( 15764 )
    if s[4] == 4:
        x.append( 24871 )
    if s[4] == 5:
        x.append( 74823 )
    if s[4] == 6:
        x.append( 134562 )
    if s[4] == 7:
        x.append( 174408 )
    if s[4] == 8:
        x.append( 354230 )
    if s[4] == 9:
        x.append( 820970 )

from numpy import *
from matplotlib.pyplot import *

gcf().close(1) # closes figure
errorbar((x),(y),yerr=yerr,marker='o')
#gca().set_yscale('log'); gca().set_xscale('log');
title("Peak Energy: "+str(s[0])); xlabel('T delay (s)'); ylabel('log(Activity)')
show()