Wednesday, April 20, 2016

Full Python Algorithm (Stripped Down from Initial NASA Code)

#Important to note that CSV_Files is a manually created directory with an arbitrary number of CSV #files
# -*- 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