# -*- coding: utf-8 -*-
# !/usr/bin/python
"""
Created on Thu Mar 17 19:08:58 2016
@author: Mike Englert
"""
import os
import numpy as np
import matplotlib.pyplot as plt
import glob
q=0
int1=[]
int2=[]
vP=[]
#indir = '\Local Disk Storage\CSV_Files'
#for root, dirs, filenames in os.walk(indir):
os.chdir("\Local Disk Storage\CSV_Files")
for file in glob.glob("*.csv"):
#q++
q=q+1
print q
#filename = '.\C1mcp&ssd00000.csv'
data = np.loadtxt(file, unpack=True, delimiter = ',', skiprows =5)
preTime = data[0]
preAmp = data[1]
#print preAmp
#Pre-Filter
preTime -= preTime[0]
preAmp2 = -preAmp
#print preAmp2
#Pass 1
iStart = 16
postAmp=preAmp2.copy()
k1 = 0.25
kn = 0.75
basePoints = 16
peakPoints = 4
for i in range (0, iStart):
v = preAmp2[i]
if i is 0:
baselineSum = mav = v
#print "baseline mav= ", v
else:
if (i < basePoints):
baselineSum += v
mav = (k1*v) + (kn*mav)
#print "mav1: ", mav
#print "baselineSum = ", baselineSum
baseline = baselineSum/basePoints
#print "baseline = ", baseline
for i in range(iStart, preAmp2.size):
mav = (k1*preAmp2[i] + (kn * mav))
#print "mav = ", mav
postAmp[i] = v = mav - baseline
if i is iStart:
vPeaks = [v for j in range(peakPoints)]
#print "VPeaks 1: ", vPeaks
iPeaks = [i for j in range(peakPoints)]
#print "IPeaks 1: ", iPeaks
else:
for j in range(peakPoints):
if v > vPeaks[j]:
vPeaks = vPeaks[:j] + [v] + vPeaks[j+1:]
#print "VPeaks 1: ", vPeaks
iPeaks = iPeaks[:j] + [i] + iPeaks[j+1:]
#print "IPeaks 1: ", iPeaks
break
#print "postAmp = ", postAmp
vPeak = sum(vPeaks)/len(vPeaks)
#print len(postAmp)
#print "Pulse Peak = ", vPeak
iPeak = round(sum(iPeaks)/len(iPeaks))
#print "Pulse Peak Index = ", iPeak
#Pass 2
part1 = 24
part2 = 270
iCFD = iStart
cfdThresh = 0.25
vCFD = vPeak * cfdThresh
integral1 = 0.0
integral2 = 0.0
for i in range(iStart, int(iPeak)):
v = postAmp[i]
if v < vCFD:
if iCFD < i:
v -= postAmp[iCFD]
integral1 += v
integral2 += v
print "integral1 = ", integral1
print "integral2 = ", integral2
#integral1=np.array(integral1)
#integral2=np.array(integral2)
iCFD += 1
else:
if i <iCFD + part1:
integral1 +=v
#integral1=np.array(integral1)
if i <iCFD + part2:
integral2 +=v
#integral2=np.array(integral2)
iPart1 = iCFD + part1
iPart2 = iCFD + part2
for i in range(int(iPeak), max(iPart1, iPart2)):
v = postAmp[i]
if i < iPart1:
integral1 +=v
#integral1=np.array(integral1)
if i < iPart2:
integral2 +=v
#integral2=np.array(integral2)
int1.append(integral1)
int2.append(integral2)
vP.append(vPeaks)
i1 = np.array(int1)
i2 = np.array(int2)
vP2 = np.array(vP)
plt.figure(1)
plt.title('Neutron PSD')
plt.xlabel('Peak Amplitude')
plt.ylabel('Integral Ratio (Long/Short)')
plt.plot(vP2, i2/i1, '.', markersize=1)
plt.figure(2)
plt.title('Integral Ratio')
plt.xlabel('Short Integral')
plt.ylabel('Long Integral')
plt.plot(i1, i2, '.', markersize=1)
plt.figure(3)
plt.title("Moving Average Smoothing Function")
plt.xlabel('Time (seconds)')
plt.ylabel('Ampitude')
plt.axis([0e0, 1.2e-06, -0.01, 0.07])
plt.plot(preTime, preAmp2, '.r', label = 'pre-Filtered Data')
plt.plot(preTime, postAmp+baseline, label = 'Smoothing Function')
plt.legend(loc = 'upper right')
No comments:
Post a Comment