Source code :: SpectrometerControl
[Return]
[Download]#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Program for controling of the Ocean optics spectrometers and extractin the most important properties from spectra
"""
from numpy import *
import matplotlib
matplotlib.use('agg')
from matplotlib.pyplot import *
import os
from time import mktime, strptime,sleep,asctime, localtime,struct_time
import re
from scipy.optimize import leastsq
from scipy.interpolate import interpolate
from scipy import stats
import scipy.stats
from scipy.special import erf
import ConfigParser
import zipfile
import socket
import shutil
import locale
import subprocess
from matplotlib.ticker import NullFormatter,ScalarFormatter
set_printoptions(precision=4,linewidth=400)
rcParams['backend'] = 'Agg'
def MAD(x):
"""
Median absolute deviation
"""
return median(abs(x-median(x)))*1.4826
def MovingAverage(data, subset_size):
"""
Simple moving average
"""
if subset_size < 1 and len(data) < subset_size:
raise ValueError('subset_size must be > 1 and less then len(data)')
mvavg = zeros(size(data))
for i in range(subset_size):
mvavg[subset_size/2:-subset_size/2]+= data[i: (+i-subset_size)]
mvavg/= float(subset_size)
mvavg[:subset_size/2] = mean(data[:subset_size/2+1])
mvavg[-subset_size/2:] = mean(data[-subset_size/2-1:])
return mvavg
def UniversalLineShape( p,x):
"""
The most general line shape. xc is position, A amplitude, w width, a skewness, n voitig profile
Line profile is integrated over single pixels
"""
x = linspace(x[0],x[-1],size(x)*10) #increase resolution
xc =p[2]
A = p[0]
w = abs(p[1])
a = p[3]
n = p[4]
# "skew voitig profile"
#http://en.wikipedia.org/wiki/Skew_normal_distribution
y = abs(A)*(1/sqrt(2*pi*w**2)*exp(-((x-xc)/w)**2/2)*(1+erf(a*(x-xc)/sqrt(2)/w))*(1-n)+ n*2/pi*w/(4*(x-xc)**2+w**2))
y = reshape(y, (size(x)/10,10) )
y = sum(y,1)/10 #integrate over pixel
return y
def lineCharacteristics(p):
"""
Calculate real properties of the line shape from function UniversalLineShape.
pc[0] - surface under line, pc[1] - FWHM, pc[2] - center of mass, pc[3] - skewness, pc[4] - voitig profile
"""
pc = zeros(5)
a = p[3]
n = p[4]
p[1] = abs(p[1])
d = a/sqrt(1+a**2)
pc[0] = p[0]
pc[2] = p[2]+(1-n)*p[1]*d*sqrt(2/pi)
f_L = p[1]/2*n
f_G = (1-n)*p[1]*sqrt(1-2*d**2/pi)*2*sqrt(2*log(2))
pc[1] = 0.5*f_L+sqrt(0.22*f_L**2+f_G**2)#http://en.wikipedia.org/wiki/Voigt_profile
pc[3] = (1-n)*(4-pi)/2*(d*sqrt(2/pi))**3/(1-2*d**2/pi)**(1.5)
pc[4] = n
return pc
class Tokamak:
"Tokamak properties and methods container class "
def __init__(self, name):
self.name = name.upper()
self.loadConfig()
self.shotNumber = -1
self.uploadServerIp = '0.0.0.0'
self.elements = list()
for name in self.ElementsList:
self.elements.append([loadtxt('SpectralLines/'+name+'.txt'),name])
def loadConfig(self):
"""
Load tokamak properties from config file
"""
if not os.path.isfile('tokamak_'+self.name+'.cfg'):
raise Exception('Tokamak config file does not exists.')
config = ConfigParser.RawConfigParser()
config.read('tokamak_'+self.name+'.cfg')
name = config.get('DEFAULT', 'Name')
if name != self.name:
raise Exception('WrongConfig')
self.dataAcqusitionTime = config.getfloat('Basic parameters', 'Plasma length')
self.trigger_shift = config.getfloat('Basic parameters', 'trigger shift')
self.database = config.get('Database', 'remote database folder')
self.trigger_port = config.getint('Database', 'trigger port')
self.login = config.get('Database', 'database login')
elements_str = config.get('Elements', 'ions')
self.ElementsList = elements_str.replace("', '", ' ').strip("'[]").split()
def saveConfig(self):
"""
Save tokamak properties to config file
"""
config = ConfigParser.RawConfigParser()
config.set('DEFAULT', 'Name', self.name)
config.add_section('Basic parameters')
config.set('Basic parameters', 'Plasma length', self.dataAcqusitionTime)
config.set('Basic parameters', 'trigger shift', self.trigger_shift)
config.add_section('Database')
config.set('Database', 'remote database folder', self.database)
config.set('Database', 'trigger port', self.trigger_port)
config.set('Database', 'database login', self.login)
config.add_section('Elements')
config.set('Elements', 'ions', self.ElementsList)
with open('tokamak_'+self.name +'.cfg', 'wb') as configfile:
config.write(configfile)
##GOLEM
def storeData(self, data_folder):
"""
Upload measured data and graphs to database
"""
#TODO odstranit!!!
#self.uploadServerIp = '192.168.2.125'
#self.shotNumber = 6879
if self.uploadServerIp == '0.0.0.0' or self.shotNumber == -1:
raise Exception('can not access to database')
if os.path.isfile(data_folder+'spectra.txt_LinesEvolution_.png'):
os.system('convert -resize 150 %sspectra.txt_LinesEvolution_.png %sicon.png' %((data_folder,)*2))
while(True):
try:
if os.system(('scp -r %s %s:%s' % (data_folder,self.login+'@'+self.uploadServerIp,self.database ))%(self.shotNumber)):
raise Exception('can not copy data')
if os.path.isfile(data_folder+'icon.png'):
if os.system(('scp -r %s %s:%s' % (data_folder+'icon.png',self.login+'@'+self.uploadServerIp,self.database+'icon.png'))%(self.shotNumber)):
raise Exception('can not copy icon')
break
except Exception as inst:
print 'problem with database access:', inst
sleep(5)
def actualShotNumber(self):
"""
Return number of the actual shot
"""
return self.shotNumber
def waitOnSoftwareTrigger(self):
"""
Waits on the signal from the main server, that database is ready. Also shot number of actual shot is reclieved.
"""
def netcat(hostname, port):
s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
s.bind((hostname, port))
s.listen(1)
conn, addr = s.accept()
data = conn.recv(1024)
conn.close()
return addr,data
HOST = '' # Symbolic name meaning all available interfaces
PORT = self.trigger_port # Arbitrary non-privileged port
#self.uploadServerIp = '192.168.2.125' #TODO vrátit!!
#self.shotNumber = 7341
#return
[addr,data] = netcat(HOST, PORT)
self.shotNumber = int(data) # je to integer
self.uploadServerIp = addr[0]
def waitOnHardwareTrigger(self, filename):
"""
Detect that spectrometer has reclieved hardware trigger
"""
#filename = './loopnumber'
#detect the hardware trigger only by change of the file with the data
if not os.path.isfile(filename):
creat_time = 0
else:
creat_time = os.stat(filename).st_mtime
change_time = creat_time
while(creat_time == change_time):
if not os.path.isfile(filename):
change_time = 0
else:
change_time = os.stat(filename).st_mtime
sleep(1)
#class Tokamak:
#def __init__(self, name):
#self.name = name.upper()
#self.loadConfig()
#self.shotNumber = -1
#self.uploadServerIp = '0.0.0.0'
#def loadConfig(self):
#if not os.path.isfile('tokamak_'+self.name+'.cfg'):
#raise Exception('Tokamak config file does not exists.')
#config = ConfigParser.RawConfigParser()
#config.read('tokamak_'+self.name+'.cfg')
#name = config.get('DEFAULT', 'Name')
#if name != self.name:
#raise Exception('WrongConfig')
#self.dataAcqusitionTime = config.getfloat('Basic parameters', 'Plasma length')
#self.database = config.get('Database', 'remote database folder')
#self.trigger_port = config.getint('Database', 'trigger port')
#self.login = config.get('Database', 'database login')
#TODO sezna očekávaných prvků! v golemovi není B, i v pořadí v jakém se nají hledat
#COMPASS
#def storeData(self, data_folder):
"""
Upload measured data and graphs to database
"""
#newpath = '../'+ str(self.shotNumber)+'/'
#if os.path.exists(newpath):
#print 'file '+newpath+'exists, will be deleted'
#shutil.rmtree(newpath)
#shutil.copytree(data_folder,newpath)
#def actualShotNumber(self):
"""
Return number of the actual shot
"""
#self.shotNumber = -1
#if os.path.exists('index.php'):
#os.remove('index.php')
#if os.system("wget --post-data='"+self.login.decode('base64')+ "' " +self.database) or not os.path.exists('index.php'):
#print 'Problems with logbook access'
#return self.shotNumber
#f = open('index.php', 'r')
#search_str = "name=shot_no class=input_text_s value="
#for i,line in enumerate(f):
#if line.find(search_str) != -1:
#ind = line.find(search_str)
#self.shotNumber = int(line[ind+len(search_str): ind+4+len(search_str)])
#break
#if self.shotNumber == -1:
#raise Exception('Problems with logbook access')
#return self.shotNumber
#def waitOnSoftwareTrigger(self):
"""
Waits on the signal from the main server, that database is ready. Also shot number of actual shot is reclieved.
"""
#pass
#def waitOnHardwareTrigger(self, filename):
"""
Detect that spectrometer has reclieved hardware trigger
"""
##detect the hardware trigger only by change of the file with the data
#if not os.path.isfile(filename):
#creat_time = 0
#else:
#creat_time = os.stat(filename).st_mtime
#change_time = creat_time
#while(creat_time == change_time):
#if not os.path.isfile(filename):
#change_time = 0
#else:
#change_time = os.stat(filename).st_mtime
#sleep(1)
class Spectrometer():
"Spectrometer properties and methods container class "
def __init__(self, SerialNumber,tokamak):
self.SerialNumber = SerialNumber
self.Tokamak = tokamak
#def fittingFunction( p,x): # ordinary gaussian
#if len(p) != 3:
#raise Exception('WrongNumberOfParameters')
#print 'spec fit fun'
#p2 = empty(5)
#p2[:3] = p
#p2[3:] = 0
#shape = UniversalLineShape( p2,x)
#return shape
#def lineCharacteristics(p): # do nothing
#return p
def elementIdentification(self,wavelength,sigma):
"""
Find the closest line to the "wavelength" from lines in tokamak database. If any line was found, amed "unknown" is returned.
"""
dispersion = (self.wavelength_range[1]-self.wavelength_range[0])/self.resolution
presicion = min(sqrt(self.wavelength_presicion**2+sigma**2),dispersion*3)
close_line_list = list()
closest_line = inf
element_name = ""
for i in self.Tokamak.elements:
index = abs(i[0]-wavelength) < presicion
if any(index):
close_line_list.append([i[0][index], i[1]])
print '\t\t\t',i[1],i[0][index], (i[0]-wavelength)[index]
close_line_i = argmin( abs(i[0]-wavelength))
close_line = i[0][close_line_i]
if abs(closest_line-wavelength)>abs(close_line-wavelength) and abs(close_line-wavelength) < presicion:
closest_line = close_line
element_name = i[1]
if isnan(closest_line):
print('\t\t\tclosest line: %4.3f , difference: %2.3f' % (wavelength, sigma))
if closest_line == inf:
closest_line = wavelength
element_name = 'unknown'
return element_name,closest_line,close_line_list
def convertCountToPhotons(self,wavelength,intensity, make_plot = False):
"""
Estimate the number of photons necessery to achieve number of counts in "intensity" on the "wavelength".
This estimation is not based on real measurement of effectivity, but only as multiplying of effectivity of CCD,
gratting and optical fiber.
"""
f1 = interpolate.interp1d(self.sensor_efficiency[:,0],self.sensor_efficiency[:,1], bounds_error=True, fill_value=0.3)
f2 = interpolate.interp1d(self.gratting_efficiency[:,0],self.gratting_efficiency[:,1], bounds_error=True, fill_value=0.1)
fiber_effectivity = 10**(-self.fiber_attenuation[:,1]/10*self.fiber_length)
f3 = interpolate.interp1d(self.fiber_attenuation[:,0],fiber_effectivity, bounds_error=True, fill_value=0.1)
#I dont't know sensitivity of CCD in UV !!
if make_plot:
plot(wavelength,self.photonsPerCount/(f1(wavelength)*f2(wavelength*f3(wavelength))))
xlabel('wavelength [nm]')
ylabel('photons/count')
show()
return intensity/(f1(wavelength)*f2(wavelength)*f3(wavelength))*self.photonsPerCount
#def estimateNoiseInCounts(self,wavelength, intensity):
#return ones(len(intensity))*self.readoutNoiseRMS
#TODO ještě přetížit ty funkce na extrakci
#TODO rovnou při extrakci určit readoutNoiseRMS
def peakFitting(self,x0,ix0,wavelength,intensity, plot_fit = False):
"""
Identify the peak around the middle in index ix0 (wavelength x0). The peak is fitted by NLLSQ and centre of mass of line is compared with lines from database. Option plot_fit will show the plot with data, their best fit and spectral lines near of this line.
"""
yn = intensity[ix0]
ixr = ix0+1
n = len(wavelength)
while ixr < n and yn > intensity[ixr]-sqrt(2):
yn = intensity[ixr]
ixr+=1
yn = intensity[ix0]
ixl = ix0
while ixl > 0 and yn > intensity[ixl-1]-sqrt(2):
yn = intensity[ixl-1]
ixl-=1
if ixr-ixl-3 <= 0: #probably false peak, too much narrow
ixr = min(n, ixr+3)
ixl = max(0, ixl-3)
a = intensity[ix0]*sqrt(2*pi*mean(self.inst_width[:,1])**2)
p0 = (a,mean(self.inst_width[:,1]),x0)
fun = lambda p,x: intensity[ixl:ixr]-self.fittingFunction(p,x)
popt,pcov,info,mesg,success = leastsq(fun, p0, args = wavelength[ixl:ixr], full_output=1)
if any(pcov == None):
pcov = ones((3,3))*nan
# calculate real paramethers of this line
popt_corr = self.lineCharacteristics(popt)
chisq=sum(info["fvec"]*info["fvec"])
doF=ixr-ixl-3
print('\t\twavelength %4.3f +/- %1.3f; chi^2/doF: %4.2f'%(popt_corr[2],sqrt(pcov[2,2]),chisq/max(doF,1)))
if chisq > scipy.stats.chi2.isf(0.1,doF): # estimate uncorrect fit
pcov*=chisq/doF
# find closest lines
line_name, line ,close_line_list= self.elementIdentification(popt_corr[2],sqrt(pcov[2,2]))
#if line == inf:
#plot_fit = True
#plot graph for easier debugging
if plot_fit:
step(wavelength[max(ixl-3,0):min(ixr+3,n)],intensity[max(ixl-3,0):min(ixr+3,n)],'r',where='mid')
step(wavelength[ixl:ixr+1]-(wavelength[ixl+1]-wavelength[ixl])/2, intensity[ixl:ixr+1],'b',where='post') # steps over the integrating area of the pixel
interv = arange(max(ixl-3,0),min(ixr+3,n))
errorbar(wavelength[interv],intensity[interv], yerr=ones(len(interv)),fmt='r.') #intensity with errorbars
errorbar(wavelength[ixl:ixr], intensity[ixl:ixr], yerr=ones(ixr-ixl),fmt='b.') #intensity used for fitting with errorbars
interv = linspace(wavelength[max(ixl-3,0)],wavelength[min(ixr+2,n-1)],100*(ixr-ixl))
plot(interv,self.fittingFunction(popt,interv), '--')
axhline(y=4, ls= '--',label = 'threshold')
print('wavelength %4.3f +/- %1.3f' % (popt_corr[2], sqrt(pcov[2,2])))
axvline(x = popt_corr[2], ls = '-.',c = 'g', label = "Center of mass")
colour = ( 'g','r','c', 'm', 'y', 'k' )
for i, element in enumerate(close_line_list):
for j,line in enumerate(element[0]):
axvline(x = line , label = element[1], c = colour[(i%len(colour))])
axis('tight')
xlabel('$\lambda$ [nm]')
ylabel('SNR[-]')
legend(loc = 'best')
title('wavelength %4.3f +/- %1.3f; chi^2/doF: %4.2f'%(popt_corr[2],sqrt(pcov[2,2]),chisq/doF))
show()
close()
d = diff(wavelength)[max(0,ix0-1)]
ixl = max(int(ix0-popt_corr[1]/d*3), ixl)
ixr = min(int(ix0+popt_corr[1]/d*3), ixr)
return popt_corr,sqrt(diag(pcov)),ixl,ixr,line_name, line
def identifyLines(self,wavelength, intensity,readoutNoiseRMS, debug = True):
"""
Detect line in the spectra given by field "intensity". Detection is based on chi2 statistics criterium. The detection threshold is set so that false detect of line in the spectrum with purely gaussian noise with sigma equal to readoutNoiseRMS is 1%. This method is more sensitive to weak lines than the naive detection of pixels higher than some threshold.
Lines are detected in infinite cycle until any line is detected. The peak is fitted by function "peakFitting" and at the and of the cycle the peaks are removed.
There are two ways how to remove peak. The peak can be substracted from the spectrum or the peak can be replaced by gaussian noise. Substrating can be used only of the instumental function is measured very precisely and therefore chi2 of the fit is close to 1. Because this is not fulfiled, tzhe second way is now used.
"""
print '\t identifing lines'
spectraParametres = list()
line_fit = list()
length = len(intensity)
derivation = zeros(length)
mean_peak_FWHM = 6 #[px] #TODO je to nesystematické, dát do vlastností spektrometru
lim = scipy.stats.chi2.isf(0.01/length,mean_peak_FWHM/2)*2 #using only positive part (50% of data), misscorrect detect probability 1%
#intensity-= median(intensity[self.black_pixels]) # dávat to sem?
#readoutNoiseRMS = MAD(intensity)
intensity/=readoutNoiseRMS
step = 0
while True:
step+= 1
if step >200:
print('infinitely cycle')
break
cumSumSpectra = cumsum((intensity*(intensity>0))**2)*2 - arange(length) #using only positive part (50% of data)
derivation[2:length-2] = cumSumSpectra[4:]-cumSumSpectra[:length-4]+mean_peak_FWHM
if all(derivation < lim):
break
#find maximum on the neighbourhood
ix0 = argmax(derivation)
ix0 = argmax(intensity[max(0,ix0-mean_peak_FWHM):min(ix0+mean_peak_FWHM, length)])+max(0,ix0-mean_peak_FWHM)
x0 = wavelength[ix0]
##if ix0 == self.resolution-1:
#plot(derivation, label= 'derivation')
##plot((0,length), (lim,lim), label = 'threshold')
#axhline( y = lim, c = 'r', ls = '--',label = 'threshold')
#axvline( x = ix0, c = 'g',label = 'max')
#axvline( x = argmax(derivation), c = 'b', ls = '--',label = 'max deriv')
##yscale('log', nonposy='clip')
#axis('tight')
#legend()
#show()
#plot(intensity, label= 'intensity')
#axvline( x = ix0, c = 'g',label = 'max')
#axvline( x = argmax(derivation), c = 'b', ls = '--',label = 'max deriv')
#axis('tight')
#legend()
#show()
if debug:
plot(wavelength, cumSumSpectra, label = 'culumative sum of squars')
axvline( x = x0, c = 'r', ls = '--')
title('cumulative sum of squars of data')
show()
popt,sig,ixl,ixr,line_name,line = self.peakFitting(x0,ix0,wavelength,intensity,debug)
#renormalize
popt[0] = self.convertCountToPhotons(x0,readoutNoiseRMS*popt[0] )
sig[0] = self.convertCountToPhotons(x0,readoutNoiseRMS*sig[0] )
sig*=stats.t.isf(0.16,ixr- ixl-3) #"1 sigma" probability (68%)
# wavelength, name, area,error,FWHM,error, Delta lambda, error
line_fit.append((line,line_name, popt[0],sig[0], popt[1],sig[1],(popt[2]-line ), sig[2]))
savetxt('current_peak.txt',(wavelength[ixl:ixr],intensity[ixl:ixr]))
# remove the fitted peak (in ideal case the peak should be substracted
#plot(wavelength[ixl+1:ixr-1], intensity[ixl+1:ixr-1],'o')
#plot(wavelength[ixl+1-5:ixr-1+5], intensity[ixl+1-5:ixr+5-1],'.r')
#show()
peakInterv = range(max(ixl+1,0),min(ixr, self.resolution-1)-1)
intensity[peakInterv] = random.randn(len(peakInterv))
return line_fit
class OceanOptics_Spec(Spectrometer):
"child of class Spectrometer with the Ocean Optics specific functions"
def __init__(self, SerialNumber,tokamak):
Spectrometer.__init__(self, SerialNumber,tokamak)
#self.SerialNumber
#self.Tokamak
#self.calibrationPolynomeCorrection = zeros((4,1))
#self.wavelength_presicion =0.1#nm
#self.black_pixels = nan
#self.resolution = 3648
#self.min_integ_time = 3.8
#self.detector = 'TCD1304AP'
#self.gratting_efficiency = array(((0,1),(1200,1)))
#self.sensor_efficiency = loadtxt('./HR4000/citlivost')
#self.photonsPerCount = 60
#self.hotPixels = None
#self.inst_width = array(((0,0),(1200,0)))
#self.inst_skewness = array(((0,0),(1200,0)))
#self.voitig_parameter = 0
self.NonLinCoeff = zeros(7) # implementovat to vůbec??
self.NonLinCoeff[0] = 1
database = os.listdir('SpectralLines')
#load everything from config file
self.loadConfig()
def start(self):
"""
It will start the java constrol program of the spectrometer.
"""
def isThisRunning( process_name ):
ps = subprocess.Popen("ps -eaf", shell=True, stdout=subprocess.PIPE)
output = ps.stdout.read()
output = output.find(process_name)
ps.stdout.close()
ps.wait()
if output == -1:
return False
else:
return True
java_library = '/opt/OmniDriver-/OOI_HOME/HighResTiming.jar:/opt/OmniDriver-/OOI_HOME/OmniDriver.jar'
if isThisRunning('HighSpeedAcquisitionSample'): #TODO dodělat
print 'spectrometer is already running'
return
try:
os.remove('nohup_java.out')
os.remove('highspeedacquisitionsample/HighSpeedAcquisitionSample.class')
except:
print 'file can not be removed'
os.system('javac -classpath '+java_library+' HighSpeedAcquisitionSample.java')
shutil.copyfile('HighSpeedAcquisitionSample.class', 'highspeedacquisitionsample/HighSpeedAcquisitionSample.class')
# need admin rights to start with so high priority
run_options = (self.SerialNumber,int(self.Tokamak.dataAcqusitionTime*1000), int(1000*self.integTime))
os.system('nohup chrt -f 99 java -Djava.library.path=. -classpath '+java_library+':. highspeedacquisitionsample.HighSpeedAcquisitionSample %s %i %i & > nohup_java.out' %run_options)
#os.system('tail -f nohup.out &')
def wavelengthCorrection(self, old_wavelength):
"""
Apply the wavelength correction polynome
"""
new_wavelegth = copy(old_wavelength)
for i,a in enumerate(self.calibrationPolynomeCorrection):
new_wavelegth+= a*old_wavelength**i
return new_wavelegth
def loadConfig(self):
"""
Load the configuration file of the spectrometer
"""
def parseStringOfField(s):
s = s.splitlines()
field = list()
for line in s:
row = (line.strip('[]')).split()
field.append(list())
for col in row:
field[-1].append(float(col))
return array(field)
if not os.path.isfile(self.SerialNumber+'.cfg'):
raise Exception('Spectrometer config file does not exists.')
config = ConfigParser.RawConfigParser()
config.read(self.SerialNumber+'.cfg')
SerialNumber = config.get('Basic parameters', 'Serial Number')
if SerialNumber != self.SerialNumber:
raise Exception('WrongSerialNumber')
self.detector = config.get('Basic parameters', 'Detector')
self.GratingTyp = config.get('Basic parameters', 'Optical grating')
self.wavelength_presicion = config.getfloat('Basic parameters', 'Wavelength presicion [nm]')
self.min_integ_time = config.getfloat('Basic parameters', 'Minimum integration time [ms]')
tmpstr = config.get('Basic parameters', 'Opticaly black pixels')
self.black_pixels = int_(parseStringOfField(tmpstr))
self.resolution = config.getint('Basic parameters', 'Resolution [px]' )
self.photonsPerCount = config.getfloat('Basic parameters', 'min photons/count')
tmpstr = config.get('Efficiencies', 'Gratting efficiency')
self.gratting_efficiency = parseStringOfField(tmpstr)
tmpstr = config.get('Efficiencies', 'Sensor efficiency') # pro HR4000 to velmi výrazně nesedí s údaji na webu
self.sensor_efficiency = parseStringOfField(tmpstr)
tmpstr = config.get('Efficiencies', 'optical fibers attenuation (dB/m)')
self.fiber_attenuation = parseStringOfField(tmpstr)
self.fiber_length = config.getfloat('Efficiencies', 'fiber length')
tmpstr = config.get('Other properties', 'Calibration polynome correction')
self.calibrationPolynomeCorrection = parseStringOfField(tmpstr)[0]
tmpstr = config.get('Other properties', 'Linearity polynome correction')
self.linearityPolynomeCorrection = parseStringOfField(tmpstr)[0]
tmpstr = config.get('Other properties', 'Hot pixels')
self.hotPixels = parseStringOfField(tmpstr)[0]
tmpstr = config.get('Other properties', 'wavelength range')
self.wavelength_range = parseStringOfField(tmpstr)[0]
tmpstr = config.get('Instrumental broadening function', 'width(wavelength)')
self.inst_width = parseStringOfField(tmpstr)
self.F_inst_width = interpolate.interp1d(self.inst_width[:,0],self.inst_width[:,1],bounds_error=False,fill_value=0)
tmpstr = config.get('Instrumental broadening function', 'skewness(wavelength)')
self.inst_skewness = parseStringOfField(tmpstr)
self.F_inst_skewness = interpolate.interp1d(self.inst_skewness[:,0],self.inst_skewness[:,1],bounds_error=False,fill_value=0)
self.voitig_parameter = config.getfloat('Instrumental broadening function', 'voitig parameter')
print 'Loaded data of ', SerialNumber
def saveConfig(self):
"""
Save the configuration file of the spectrometer
"""
config = ConfigParser.RawConfigParser()
config.add_section('Basic parameters')
config.set('Basic parameters', 'Type', 'Ocean Optics Spectrometer')
config.set('Basic parameters', 'Serial Number', self.SerialNumber )
config.set('Basic parameters', 'Detector', self.detector )
config.set('Basic parameters', 'Optical grating', self.GratingTyp )
config.set('Basic parameters', 'Date', asctime())
config.set('Basic parameters', 'Wavelength presicion [nm]', self.wavelength_presicion)
config.set('Basic parameters', 'Opticaly black pixels', self.black_pixels )
config.set('Basic parameters', 'Resolution [px]', self.resolution )
config.set('Basic parameters', 'Minimum integration time [ms]', self.min_integ_time)
config.set('Basic parameters', 'min photons/count', self.photonsPerCount)
config.add_section('Efficiencies')
config.set('Efficiencies', 'Gratting efficiency', self.gratting_efficiency )
config.set('Efficiencies', 'Sensor efficiency', self.sensor_efficiency)
config.set('Efficiencies', 'Optical fibers attenuation (dB/m)', self.fiber_attenuation)
config.set('Efficiencies', 'fiber length', self.fiber_length)
config.add_section('Other properties')
config.set('Other properties', 'Calibration polynome correction', self.calibrationPolynomeCorrection)
config.set('Other properties', 'linearity polynome correction' , self.linearityPolynomeCorrection)
config.set('Other properties', 'Hot pixels', self.hotPixels)
config.set('Other properties', 'wavelength range',self.wavelength_range )
config.add_section('Instrumental broadening function')
config.set('Instrumental broadening function', 'width(wavelength)', self.inst_width ) # šířka kalibračních čar v několika bodech
config.set('Instrumental broadening function', 'skewness(wavelength)', self.inst_skewness ) # šikmost kalibračních čar v několika bodech
config.set('Instrumental broadening function', 'voitig parameter', self.voitig_parameter ) # negausovkost základny
with open(self.SerialNumber+'.cfg', 'wb') as configfile:
config.write(configfile)
def checkCalibration(self): # vykresí to graf a určí chi^2 = suma(delta labda ku sigma)^2
#zapolá to makeCalibration
pass
def makeCalibration(self, dryRun = True):
pass
#načte list kalibračních čar
#100x spustí sběr data a zprůměruje
#projde kalibrační čáry a najde nejbližší peak
#aktualizuje (pro dryRun = false) to self.inst_width, self.inst_skewness, calibrationPolynomeCorrection,
# vykreslí to současné a nové hodnoty (i errorbary)
#TODO nastavit width(wavelength),skewness, voitig parameter, Calibration polynome correction
def calcParamsForFitting(self,p):
"""
The lines are not fitted by full set of paramathes of the general peak function UniversalShape. Only a few paramethers are used, other are calculated here.
"""
p_full = empty(5)
p_full[:3] = p[:3]
#print 'int_w ', self.F_inst_width(p[2]), self.F_inst_skewness(p[2]), self.voitig_parameter
#p_full[1] = p[1]
p_full[1] = sqrt(p[1]**2+self.F_inst_width(p[2])**2) #add instrumental broadening
p_full[3] = self.F_inst_skewness(p[2])
#p_full[3] = 10
p_full[4] = self.voitig_parameter
#p_full[4] = 0.00
return p_full
def fittingFunction(self,p,x): # využívat informace o znalosti instrumentálního rozšíření
p_full = self.calcParamsForFitting(p)
return UniversalLineShape( p_full,x)
def lineCharacteristics(self,p):
"""
The lines are not fitted by full set of paramathes of the general peak function UniversalShape.
Therefore the real properties of the peak are calcutelad here from the limited set of the paramethers.
"""
p_full = self.calcParamsForFitting(p)
p_corrected = lineCharacteristics(p_full)
p_corrected[1] = sqrt(p_corrected[1] **2-self.F_inst_width(p_corrected[2])**2)
return lineCharacteristics(p_full)
def loadData(self,path,typ,name = 'spectra'):
"""
Load data from file and create object "shot" containig all necessary informations andou the tokamak shot. It supports many different types of files.
Java_Data_file - file experted from java control program HighSpeedAcquisitionSample.java
Python_Data_file - file exported from this python program
SpectraSuite_xml - file exported from OceanOptics SpectralSuite
SpectraSuite_txt - acsii file exported from OceanOptics SpectralSuite
SpectraSuite_HighSpeed - ascii file exported from OceanOptics SpectralSuite after high speed acquisition
"""
shot = Shot(self.Tokamak.shotNumber, self)
if typ == 'Java_Data_file': # acqurired by HighSpeedAcquisitionSample.java
#extract data about spectra
shot.DataFile = path+name+'.txt'
if not os.path.isfile(shot.DataFile):
raise Exception('Data file '+shot.DataFile+' does not exists.')
f = open(shot.DataFile, 'r')
for i,line in enumerate(f):
if line[:14] == 'Date and time:':
shot.time = mktime(strptime(line[15:-1], "%Y.%m.%d %H:%M:%S"))
if line[:18] == 'Number of spectra:':
shot.n_spectra = int(line[20:-1])
if line[:25] == 'Number pixels in spectra:':
shot.n_pixels = int(line[27:-1])
if line[:22] == 'Integration time [us]:':
shot.integ_time = int(line[24:-1])/1000.0
if line[:22] == 'Board temperature [C]:':
shot.temperature = float(line[24:-1])
if line[:17] == 'Spectrometer S/N:':
SN = line[18:-1]
if line[:18] == 'Exact time stamps:':
break
if SN != self.SerialNumber:
raise Exception('wrong serial number, '+SN+' != '+self.SerialNumber)
shot.time_stamps = zeros(shot.n_spectra)
for i,line in enumerate(f):
if i >= shot.n_spectra:
break
m = re.search('\s[0-9]+\.[0-9E]+',line)
s = m.group(0)
shot.time_stamps[i] = float(s)
shot.time_stamps/= 1e6 #convert ns -> ms
shot.time_stamps+= -shot.time_stamps[0]+self.Tokamak.trigger_shift #shift ->0
data = genfromtxt(f)
f.close()
shot.wavelength = data[:shot.n_pixels]
shot.readoutNoiseRMS = 0
for i in range(shot.n_spectra):
intensity = data[(i+1)*shot.n_pixels: (i+2)*shot.n_pixels]
#intensity -= median(intensity[self.black_pixels])
shot.spectra.append(Spectrum(shot,shot.wavelength,intensity,shot.time_stamps[i] ))
shot.readoutNoiseRMS += std(intensity[self.black_pixels])**2
shot.readoutNoiseRMS = sqrt(shot.readoutNoiseRMS/shot.n_spectra) # Zkontrolovat stabilitu odhadu
#extract data about background and dark current
if os.path.isfile(path+'background.txt'):
imax = 0
f = open(path+'background.txt', 'r')
for i,line in enumerate(f):
if line[:-1].replace('.','').isalnum():
imax = i
break
if line[0:22] == 'Integration time [us]:':
integ_time_background = int(line[24:-1])/1000
f.close()
shot.darkNoise = genfromtxt(path+'background.txt', skip_header=imax)
shot.darkNoise -= median(shot.darkNoise[self.black_pixels])
shot.darkNoise*=shot.integ_time/float(integ_time_background)
#extract data about readout patterns
if os.path.isfile(path+'readoutpatterns.txt'):
imax = 0
f = open(path+'readoutpatterns.txt', 'r')
for i,line in enumerate(f):
if line[:-1].replace('.','').isalnum():
imax = i
break
f.close()
shot.readOutPatterns = genfromtxt(path+'readoutpatterns.txt', skip_header=imax)
shot.readOutPatterns-= (shot.darkNoise/shot.integ_time)*self.min_integ_time
else:
shot.readOutPatterns = mean(intensity[self.black_pixels])
if typ == 'Python_Data_file':
shot.DataFile = path+name+'.txt'
if not os.path.isfile(shot.DataFile):
raise Exception('Data file '+shot.DataFile+' does not exists.')
f = open(shot.DataFile, 'r')
for i,line in enumerate(f):
ind = line.find('Serial Number')
if ind != -1:
SN = line[ind+15:-1]
ind = line.find('Date and time (GMT)')
if ind != -1:
shot.time = mktime(strptime(line[ind+21:-1]))
ind = line.find('Number of spectra')
if ind != -1:
shot.n_spectra = int(line[ind+19:-1])
ind = line.find('Resolution')
if ind != -1:
shot.n_pixels = int( line[ind+16:-1])
ind = line.find('Integration time [ms]')
if ind != -1:
shot.integ_time = float(line[ind+22:-1])
ind = line.find('Board temperature [C]')
if ind != -1:
shot.temperature = float(line[ind+22:-1])
ind = line.find('Time stamps [ms]')
if ind != -1:
s = line[ind+18:-1]
shot.time_stamps = float_(s.strip('[]').split())
ind = line.find('Noise RMS')
if ind != -1:
print line[ind+20:-1]
shot.readoutNoiseRMS = float(line[ind+20:-1])
if line.find('***************') != -1:
break
if SN != self.SerialNumber:
raise Exception('wrong serial number, '+SN+' != '+self.SerialNumber)
data = genfromtxt(f)
f.close()
shot.wavelength = data[:,0]
for i in range(shot.n_spectra):
shot.spectra.append(Spectrum(shot,shot.wavelength, data[:,i+1], shot.time_stamps[i]))
shot.darkNoise = zeros(shot.n_pixels)
shot.readOutPatterns = zeros(shot.n_pixels)
if typ == 'SpectraSuite_xml': #acqurid by Data Studio
shot.temperature = nan
shot.DataFile = path+name+'.ProcSpec'
if not os.path.isfile(shot.DataFile):
raise Exception('Data file '+shot.DataFile+' does not exists.')
if os.path.isfile('./readoutpatterns.txt'):
shot.readOutPatterns = genfromtxt('./readoutpatterns.txt', skip_header=1)
shot.readOutPatterns -= median(shot.readOutPatterns[self.black_pixels])
with zipfile.ZipFile(shot.DataFile, 'r') as myzip:
shot.n_spectra = len(myzip.namelist())-2
for data in myzip.namelist():
if data != 'OOISignatures.xml' and data != 'OOIVersion.txt':
#print directory+'/'+data_file+'/'+data
t = myzip.getinfo(data).date_time
shot.time = mktime((t[0],t[1],t[2],t[3],t[4],t[5],0, 0,0))
lines = myzip.open(data)
for j,line in enumerate(lines):
if line.find('<numberOfPixels>') != -1:
ind = line.find('<numberOfPixels>')
shot.n_pixels = int(line[ind+16:-18])
if line.find('<integrationTime>') != -1:
ind = line.find('<integrationTime>')
shot.integ_time = int(line[ind+17:-19])/1000.0#convert to ms
if line.find('<spectrometerSerialNumber>') != -1:
ind = line.find('<spectrometerSerialNumber>')
SN = line[26+ind:-28]
if SN != self.SerialNumber:
raise Exception('wrong serial number, '+SN+' != '+self.SerialNumber)
lines.close()
lines = myzip.open(data)
j = 0
intensity = zeros(shot.n_pixels)
shot.wavelength = zeros(shot.n_pixels)
channelWavelengths = False
sourceSpectra = False
pixelValues = False
for i,line in enumerate(lines):
index = line.find('<double>')
if line.find('<com.oceanoptics.omnidriver.spectra.OmniSpectrum>') != -1:
sourceSpectra = True
j = 0
if line.find('</com.oceanoptics.omnidriver.spectra.OmniSpectrum>') != -1:
sourceSpectra = False
if line.find('<channelWavelengths>') != -1:
channelWavelengths = True
j = 0
if line.find('</channelWavelengths>') != -1:
channelWavelengths = False
if line.find('<pixelValues>') != -1:
pixelValues = True
j = 0
if line.find('</pixelValues>') != -1:
pixelValues = False
if index != -1:
if sourceSpectra and pixelValues :
intensity[j] = float(line[index+8:-10])
if sourceSpectra and channelWavelengths:
shot.wavelength[j] = float(line[index+8:-10])
j+=1
lines.close()
intensity -= median(intensity)
shot.spectra.append(Spectrum(shot,shot.wavelength,intensity))
shot.time_stamps = arange(shot.n_spectra)*shot.integ_time+self.Tokamak.trigger_shift
shot.readoutNoiseRMS = MAD(intensity)
if typ == 'SpectraSuite_txt': #text acqurid by Data Studio
shot.temperature = nan
shot.DataFile = path+name+'.txt'
if not os.path.isfile(shot.DataFile):
raise Exception('Data file '+shot.DataFile+' does not exists.')
if os.path.isfile('./readoutpatterns.txt'):
shot.readOutPatterns = genfromtxt('./readoutpatterns.txt', skip_header=1)
shot.readOutPatterns -= median(shot.readOutPatterns[self.black_pixels])
f = open(shot.DataFile, 'r')
line = f.readline()
if line != 'SpectraSuite Data File\n':
raise Exception('wrong file type')
while True:
line = f.readline()
if not line:
if not line: raise Exception('corrupted header of file')
if line[:5] == 'Date:':
try:
shot.time = mktime(strptime(line[6:-1], "%a %b %d %H:%M:%S %Z %Y")) # sometimes it can make mysterious errors connected with figure()
except:
shot.time = os.path.getctime(shot.DataFile)
if line[:36] == 'Number of Sampled Component Spectra:':
shot.n_spectra = int(line[37:-1])
if line[:16] == 'Number of Pixels':
shot.n_pixels = int(line[-5:-1])
if line[:24] == 'Integration Time (usec):':
if line[-11:-1] == '('+self.SerialNumber+')':
shot.integ_time = int(line[25:-11])/1000.0 #convert to ms
else:
shot.integ_time = int(line[25:-1])/1000.0 #convert to ms
if line[:14] == 'Spectrometers:' or line[:27] == 'Spectrometer Serial Number:': #TODO opravit od do
SN = line[-9:-1]
if SN != self.SerialNumber:
raise Exception('wrong serial number, '+SN+' != '+self.SerialNumber)
if line[:5] == '>>>>>':
break
if shot.n_spectra == 0: # missing information in file
shot.n_spectra = 1
replaceComma = lambda s:float(s.replace(',','.'))
data = genfromtxt(f, skip_footer = 1, converters = {0:replaceComma, 1:replaceComma })
f.close()
shot.wavelength = data[:,0]
intensity = data[:,1]
intensity -= median(intensity)
shot.spectra.append(Spectrum(shot,shot.wavelength,intensity))
shot.readoutNoiseRMS = MAD(intensity)
#shot.readoutNoiseRMS = std(intensity[self.black_pixels])
if typ == 'SpectraSuite_HighSpeed':
shot.temperature = nan
shot.DataFile = path+name+'.txt'
if not os.path.isfile(shot.DataFile):
raise Exception('Data file '+shot.DataFile+' does not exists.')
f = open(shot.DataFile, 'r')
first_line = f.readline()
if first_line == 'SpectraSuite Data File\n':
raise Exception('wrong file type')
data = loadtxt(f)
f.close()
shot.time_stamps = float_(first_line.split())+self.Tokamak.trigger_shift
shot.time = os.path.getctime(shot.DataFile) # it can be wrong time!
shot.integ_time = shot.time_stamps[-1]/(len(shot.time_stamps)-1) # only estimation!
shot.n_spectra = size(data,1)-1
shot.n_pixels = size(data,0)
shot.wavelength = data[:,0]
for i in range(shot.n_spectra):
intensity = data[:,i+1]
intensity -= median(intensity)
shot.spectra.append(Spectrum(shot,shot.wavelength,intensity,shot.time_stamps[i] ))
shot.readoutNoiseRMS = MAD(data[self.black_pixels,1:])
if max(abs(amin(shot.wavelength)-self.wavelength_range[0]), abs(amax(shot.wavelength)-self.wavelength_range[1]))>10:
raise Exception('bad wavelength range %d - %d != %d - %d' %(amin(shot.wavelength),amax(shot.wavelength),self.wavelength_range[0],self.wavelength_range[1]))
if shot.n_pixels != self.resolution:
raise Exception('wrong number of pixels, '+str(shot.n_pixels)+' != '+str(self.resolution))
return shot
##TODO ve fitovací funkci zahrout i to , že je to přenormované šumem. (prostě vypočítávat šum i pro tu teoretickou funkci)
#class CIII_HR_Spec(Spectrometer):
##TODO různé konfigy prop různá nastavení?
#def __init__(self):
#Spectrometer.__init__(self, SerialNumber,tokamak)
##TODO musí si to pamatovat všechny nastavitelné parametry
#def convertCountToPhotons(self,wavelength, intensity, make_plot = False):
##TODO tady by si to mohlo pamatovat ten šum
##TODO založit to na vzorci z výzkumáku
#if len(wavelength) != len(intensity):
#wavelength = ones(len(intensity))*mean(wavelength)
#if make_plot:
#plot(wavelength,f1(wavelength)*f2(wavelength)*self.photonsPerCount)
#xlabel('wavelength [nm]')
#ylabel('photons/count')
#show()
##TODO střeva
#return corr_intensity # šum i intenzita přepočtené na fotony
#def estimateNoiseInCounts(self,wavelength, intensity):
#if len(wavelength) != len(intensity):
#wavelength = ones(len(intensity))*mean(wavelength)
##TODO dodělat střeva
#return noise #šum v počtu countů
#def PeakFitting:
#def fittingFunction:# buď
## přidá navíc paramer základna, využije se tu maxumum dosavadních znalostí.
#def identifyLines(self,wavelength, intensity):# return table of lines, wavelngths, +other parametres
class Shot():
"""Discharge container class
used for plotting and obtaining discharge data
"""
def __init__(self, shot_num, Spectrometer):
self.Spectrometer = Spectrometer
self.shot_num=shot_num
self.spectra = list()
self.wavelength = nan # udělat korekci při načítání
self.readOutPatterns = zeros(Spectrometer.resolution) #odečíst mediaán!!!
self.darkNoise = zeros(Spectrometer.resolution) #odečíst mediaán!!!, normované na integrační dobu
self.readoutNoiseRMS = 13 # TODO co to je? chovám se k tomu jak o k poli
self.maxSNR = 0 #TODO co stím?
self.n_spectra = 0
self.time_stamps = zeros(1)
self.DataFile = ''
def plasmaShot(self):
"""
Detect if there was any plasma during the shot.
"""
plasma = zeros(self.n_spectra, dtype = 'bool')
for i in arange(self.n_spectra):
plasma[i] = self.spectra[i].plasma()
return any(plasma)
def getData(self,sensitivity_correction = False,noiseCorrected = True):
"""
Return the field containig spectra from the all acqusitions during the shot. Option sensitivity_correction set that the correction on the not constant sensitivity of the whole device at different wavelength shold by applied. Option noiseCorrected set that the readout noise patterns, dark noise and hot pixels should be removed.
"""
spectra = empty((self.Spectrometer.resolution, self.n_spectra))
for i in arange(self.n_spectra):
spectra[:,i] = self.spectra[i].getData(sensitivity_correction , noiseCorrected )
return copy(self.wavelength), spectra
def saveData(self,filename = None, noiseCorrected = True):
"""
Export data from the shot in well arranged shape. At the bigining of the file is the header and than coloumns with wavelengths and intensities.
"""
print 'saving data '#, (self.time_stamps[-1]/(len(self.time_stamps)-1))
if filename == None:
filename = self.DataFile+'_Data_.txt'
# make a header
header = ''
header += 'Spectrometer:\t' + 'Ocean Optics Spectrometer\n'
header += 'Serial Number:\t' + self.Spectrometer.SerialNumber+'\n'
header += 'Date and time (GMT):\t' + asctime(localtime(self.time))+'\n'
header += 'Number of spectra:\t'+ str(self.n_spectra)+'\n'
header += 'Resolution [px]:\t' + str(self.Spectrometer.resolution)+'\n'
header += 'Integration time [ms]:\t'+ '%2.2f' %(self.integ_time)+'\n'
header += 'Acquisition time [ms]:\t'+ '%2.2f' %((self.time_stamps[-1]-self.time_stamps[0])/(len(self.time_stamps)-1))+'\n'
header += 'Board temperature [C]:\t'+ '%2.3f' %self.temperature +'\n'
header += 'Time stamps [ms]:\t' + str(self.time_stamps)+'\n'
header += 'Noise RMS [counts]:\t' + '%2.1f' %self.readoutNoiseRMS+'\n'
header += '*'*100+'\n'
#save spectra
spectra = empty((self.Spectrometer.resolution, self.n_spectra+1))
spectra[:,0] = self.wavelength
for i in arange(self.n_spectra):
spectra[:,i+1] = self.spectra[i].getData(noiseCorrected = noiseCorrected)
f = open(filename, 'w')
f.write(header)
savetxt(f, spectra,fmt=['%1.3f']+['%4i']* self.n_spectra)
f.close()
def plot(self, file_type = 'pdf', log_scale = True ,sensitivity_correction = True, w_interval = None):
"""
Plot spectra in separate files.
Options:
file_type - file type of exported graph
log_scale - set log scale
sensitivity_correction - correction on the not constant sensitivity of the whole device at different wavelength
w_interval - interval if the plotted wavelengths
"""
print 'plotting separate graphs'
for i,spectrum in enumerate(self.spectra):
self.spectra[i].plot(sensitivity_correction, log_scale =log_scale, noiseCorrected = True, w_interval = w_interval)
if sensitivity_correction:
ylabel('photons [-]')
else:
ylabel('counts [-]')
if file_type == None:
show()
else:
savefig(self.DataFile+'_Graph_%d.'%(i)+file_type)
close()
def plot_all(self, file_type = 'pdf', log_scale = True ,sensitivity_correction = True, w_interval = None, plasma = True):
"""
Plot spectra in to one plot
Options:
file_type - file type of exported graph
log_scale - set log scale
sensitivity_correction - correction on the not constant sensitivity of the whole device at different wavelength
w_interval - interval if the plotted wavelengths
plasma - only spectra with some lines
"""
print 'plotting all together'
class MyFormatter(ScalarFormatter):
def __call__(self, x, pos=None):
if pos==0:
return ''
else: return ScalarFormatter.__call__(self, x, pos)
fig = figure(num=None, figsize=(12, 12), dpi=80, facecolor='w', edgecolor='k')
subplots_adjust(wspace=0,hspace=0)
if not self.plasmaShot():
plasma = False
n_plots = self.n_spectra
for i,spectrum in enumerate(self.spectra):
if not spectrum.plasma() and plasma:
n_plots-= 1
i_plot = 0
for i,spectrum in enumerate(self.spectra):
if (not spectrum.plasma()) and plasma:
continue
i_plot+= 1
ax = fig.add_subplot(n_plots,1,i_plot)
ax.xaxis.set_major_formatter( NullFormatter() )
ax.yaxis.set_major_formatter( MyFormatter() )
#subplot( self.n_spectra ,1, i+1)
self.spectra[i].plot(sensitivity_correction, log_scale =log_scale, noiseCorrected = True, w_interval = w_interval)
if sensitivity_correction:
ylabel('photons [-]', rotation='horizontal')
else:
ylabel('counts [-]', rotation='horizontal')
ax = fig.add_subplot(n_plots ,1,n_plots)
ax.xaxis.set_major_formatter( ScalarFormatter() )
if file_type == None:
show()
else:
savefig(self.DataFile+'_Graph_.'+file_type)
close()
def processShot(self, details = False):
"""
prepare table of identified lines and their properties for another functions
"""
LinesTable = []
def addNewLine(wavelength,name):
LinesTable.append(list((wavelength,name)))
if details:
emptyCell = list((0,nan,nan,nan,nan,nan))
else:
emptyCell = list((0,))
for i in range(self.n_spectra):
LinesTable[-1].append(emptyCell)
def setValue(wavelength,time_slice, data):
for row in LinesTable:
if row[0] == wavelength:
if details:
row[2+time_slice] = data
else:
row[2+time_slice] = data[0]
break
for i,spectrum in enumerate(self.spectra):
Table = spectrum.processSpectrum()
print 'processing '+str(i)+'-th spectrum'
for line in Table:
exist = False
for row in LinesTable:
if row[0] == line[0]:
exist = True
setValue(line[0],i,line[2:] )
break
if not exist:
addNewLine(line[0],line[1])
setValue(line[0],i,line[2:] )
ColoumnNames = list(('wavelength','line'))
if details:
OtherColoumns = list(('area','err','FWHM','err', 'Delta lambda','err'))
else:
OtherColoumns = list(('area',))
for i in range(self.n_spectra):
for name in OtherColoumns:
ColoumnNames.append(name)
return ColoumnNames,LinesTable,self.time_stamps
def plotLinesTimeEvolution(self,file_type = None, only_identified = True, n_lines = 6):
"""
Identify annd plot the most intensive lines and plot their time evolution during the shot
Options:
file_type - file type of exported graph
only_identified - extract only identified lines
n_lines - number of plotted lines
"""
print 'lines evolution plotting'
IntensitiesTable = []
ErrorsTable = []
NamesTable = []
ColoumnNames,LinesTable,time_stamps = self.processShot(True)
for line in LinesTable:
if only_identified == True and line[1] == 'unknown':
continue
NamesTable.append(line[1]+' %3.1f' %line[0]+'nm')
intensity = zeros(self.n_spectra)
errorbars= zeros(self.n_spectra)
for j, time_slice in enumerate(line[2:]):
if not self.spectra[j].plasma():
continue
intensity[j] = time_slice[0]
errorbars[j] = time_slice[1]
IntensitiesTable.append(intensity)
ErrorsTable.append(errorbars)
if len(IntensitiesTable) == 0:
return
plasma_end = self.n_spectra
for i in arange(self.n_spectra,0,-1):
if self.spectra[i-1].plasma():
plasma_end = i
break
IntensMax = amax(IntensitiesTable, axis = 1)
sort_ind = argsort(IntensMax, order=None)[::-1]
for i,index in enumerate(sort_ind):
if i > n_lines:
break
errorbar(time_stamps,IntensitiesTable[index],yerr=ErrorsTable[index],label=NamesTable[index])
legend(loc = 'best')
t_start= time_stamps[0]
t_end = time_stamps[min(plasma_end,self.n_spectra-1)]
axis([t_start,t_end,0, amax(IntensitiesTable) *1.1])
xlabel('t [ms]')
ylabel('Intensity [a.u.]')
title('The most intense lines')
if file_type == None:
show()
else:
savefig(self.DataFile+'_LinesEvolution_.'+file_type)
close()
def positionHistorogram(self):
"""
plot the hytorogram of positions of identified lines. It is useful for false line identifications.
"""
ColoumnNames,LinesTable,time_stamps = self.processShot(True)
hist_list = list()
for line in LinesTable:
for cell in line[2:]:
if not isnan(cell[4]):
hist_list.append(( line[0] ,cell[4]))
hist_list = array(hist_list)
plot(hist_list[:,0], hist_list[:,1], '.')
show()
def tableOfResultes(self,sufix = '', sort = False, details = False, only_identified = True):
"""
Export table with lines and their properties
"""
ColoumnNames,LinesTable,time_stamps = self.processShot(details)
f = open(self.DataFile+'_Resultes_'+sufix+'.csv', 'w')
for name in ColoumnNames:
f.write(name)
f.write(' ;')
wavelengths = []
for i,row in enumerate(LinesTable):
wavelengths.append(row[0])
if sort:
wavelengths.sort()
for pos in wavelengths: # complicate way to sort the lines
for row in LinesTable:
if only_identified and row[1] == 'unknown':
continue
if row[0] == pos:
f.write('\n')
for cells in row:
col = 0
if isinstance(cells,list) or isinstance(cells,tuple):
for num in cells:
f.write(locale.format("%5.3f",num)+' ;')
col += 1
else:
if isinstance(cells, double):
f.write(locale.format('%3.3f', cells)+' ;')
else:
f.write(str(cells)+' ;')
col += 1
f.closed
class Spectrum():
"""
Spectrum container class used for plotting and preparing of single spectrum
"""
def __init__(self,Shot,wavelength,intensity, time_stamp = 0):
self.Shot = Shot
self.intensity = intensity #odečíst mediaán!!!, kontrolovat jeszli to neovlivňují další funkce
self.wavelength = Shot.Spectrometer.wavelengthCorrection(wavelength)
self.time_stamp = time_stamp
self.lineFitList = []
def plasma(self):
"""
Detect if in the data are some lines
"""
if any(self.getData() > self.Shot.readoutNoiseRMS*5):
return True
return False
def plot(self,sensitivity_correction = True, log_scale = True, noiseCorrected = True, w_interval = None):
"""
Plot a single spectrum
log_scale - set log scale
sensitivity_correction - correction on the not constant sensitivity of the whole device at different wavelength
w_interval - interval if the plotted wavelengths
"""
if not w_interval:
w_interval = [amin(self.wavelength), amax(self.wavelength)]
w_interval_bool = (self.wavelength >= w_interval[0]) * (self.wavelength <= w_interval[1] )
intensity = self.getData(sensitivity_correction,noiseCorrected)*w_interval_bool
noise = self.Shot.readoutNoiseRMS
threshold = scipy.stats.norm.isf(0.05/self.Shot.Spectrometer.resolution) # 5% probability of mistake
if sensitivity_correction:
noise = self.Shot.Spectrometer.convertCountToPhotons(self.wavelength, noise)
else:
noise *= ones(self.Shot.Spectrometer.resolution)
min_threshold = amin(threshold*noise)
noise *= w_interval_bool
max_threshold = amax(threshold*noise)
plot(self.wavelength, intensity,'b', label= 'time %2.1f ms' %self.time_stamp, linewidth=0.5)
plot(self.wavelength,threshold*noise, 'g')
xlabel('$\lambda$ [nm]',fontdict={'fontsize':12})
if log_scale:
yscale('log', nonposy='clip')
axis([w_interval[0],w_interval[1],min_threshold/8, max(max_threshold,amax(intensity))])
else:
axis([w_interval[0],w_interval[1],amin(-noise) , max(max_threshold,amax(intensity))])
legend(loc='best')
def getData(self,sensitivity_correction = False,noiseCorrected = True):
"""
prepare intensities for following processing.
In the fist step nonlinearity is removed, than if noiseCorrected is allowed the readout patterns, dark noise and hot pixels are removed and finally if sensitivity_correction is allowed the different sentivity of whole device for different wavelengths is corrected.
"""
#plot(self.intensity)
#plot( self.Shot.readOutPatterns)
#plot(self.Shot.darkNoise)
#show()
intensity = copy(self.intensity)
if noiseCorrected:
intensity -= self.Shot.readOutPatterns+self.Shot.darkNoise
intensity[int_(self.Shot.Spectrometer.hotPixels)] = median(intensity)
correction = zeros(size(intensity))
for i,a in enumerate(self.Shot.Spectrometer.linearityPolynomeCorrection[::-1]):
correction*= intensity
correction += a
intensity = intensity/correction
if sensitivity_correction:
intensity = self.Shot.Spectrometer.convertCountToPhotons(self.wavelength, intensity)
return intensity
def processSpectrum(self):
#if not self.plasma():
#return list()
if len(self.lineFitList) != 0:
return self.lineFitList
corr_intensity = self.getData(False,True)
#print 'identifyLines'#, self.Shot.readoutNoiseRMS
self.lineFitList = self.Shot.Spectrometer.identifyLines(self.wavelength,corr_intensity,self.Shot.readoutNoiseRMS, False ) # return table of lines, wavelength, +other parametres
return self.lineFitList
def main():
#print 'spektrometr'
#x = linspace(0,1,500)
#p = empty(5)
#p[0] = 1
#p[1] = 0.01
#p[2] = 0.5
#p[3] = 0.2
#p[4] = 0.15
#line = UniversalLineShape( p,x)
#print lineCharacteristics(p)
#print sum(line/500)
#u = sum(x*line)/sum(line)
#print sqrt(sum((x-u)**2*line)/sum(line))*2*sqrt(2*log(2))
#print sum(x*line)/sum(line)
#s = sqrt(sum((x-u)**2*line)/sum(line))
#print (sum((x-u)**3*line)/sum(line))/s**3
#plot(x,line)
#show()
#from SpectrometerControl import *
#COMPASS = Tokamak('COMPASS')
#COMPASS.saveConfig()
##exit()
##Spectrometer('gfdg',GOLEM)
##HR4000 = OceanOptics_Spec('HR4C1947',GOLEM)
##HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
#GOLEM = Tokamak('GOLEM')
#HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
#HR2000 = OceanOptics_Spec('HR+C1103',COMPASS)
#shot = HR2000.loadData('./data/', 'Java_Data_file', 'spectra')
#shot.saveData(noiseCorrected = True)
#shot.plot_all( None, log_scale = True, sensitivity_correction = True)
##for spec in shot.spectra:
##shot.spectra[0].intensity+= spec.intensity
##shot.readoutNoiseRMS*= sqrt(shot.n_spectra)
##shot.n_spectra = 1
##shot.tableOfResultes(details = True, sort = False, only_identified = True)
##shot.plotLinesTimeEvolution(file_type = None, only_identified = True, n_lines = 50)
#shot.positionHistorogram()
#exit()
#shot = HR2000.loadData('./data_VW/','simple','spectra')
#shot.tableOfResultes(details = False, sort = True, only_identified = True)
#shot.plot_all(None)
#shot.plotLinesTimeEvolution('png',only_identified = True)
#shot.spectra[2].plot(sensitivity_correction = True, log_scale = False, noiseCorrected = True, w_interval = None)
#show()
#hist(shot.spectra[1].intensity,70, range = (-40, 100))
#show()
#exit()
#shot.getData()
#shot.plot( None, log_scale = False, sensitivity_correction = True )
#plot(shot.darkNoise)
#show()
#shot = HR2000.loadData('./data/','SpectraSuite_HighSpeed','shotd')
#GOLEM = Tokamak('GOLEM')
#HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
#shot = HR2000.loadData('./','Python_Data_file','spectra.txt_Data_')
#shot.plasmaShot()
#shot.plot_all(None, log_scale = False, sensitivity_correction = True, plasma = True)
##shot.tableOfResultes(details = True, sort = True, only_identified = False)
#shot.plotLinesTimeEvolution('png',only_identified = True)
##show()
#colour = ( 'g','r','c', 'm', 'y', 'k' )
#shot.spectra[1].plot( log_scale = False, sensitivity_correction = True)
#for i, element in enumerate(GOLEM.elements):
#for j,line in enumerate(element[0]):
#axvline(x = line , c = colour[(i%len(colour))])
#show()
#for i, element in enumerate(GOLEM.elements):
#shot.spectra[1].plot( log_scale = False, sensitivity_correction = True)
#for j,line in enumerate(element[0]):
#axvline(x = line , c = colour[(i%len(colour))])
#title(element[1])
#show()
#exit()
#COMPASS = Tokamak('COMPASS')
#HR2000 = OceanOptics_Spec('HR+C0732',COMPASS)
#shot = HR2000.loadData('./data_nova/02.12.2011/','SpectraSuite_xml','shot000')
#for i in range(1,shot.n_spectra):
#shot.spectra[0].intensity += shot.spectra[i].intensity
#shot.spectra[0].intensity/= shot.n_spectra
#shot.readoutNoiseRMS/= sqrt(shot.n_spectra)
#shot2 = HR2000.loadData('./data_nova/02.12.2011/','SpectraSuite_xml','shot001')
#shot2.spectra[0].intensity = shot.spectra[0].intensity
#shot2.plot_all( None, log_scale = False, sensitivity_correction = True)
#shot2.tableOfResultes(details = True, sort = True, only_identified = False)
#shot.saveData(noiseCorrected = True)
#exit()
#shot.spectra[0].intensity = random.randn(HR2000.resolution)*shot.readoutNoiseRMS
#plot(shot.spectra[0].intensity)
#show()
#shot.tableOfResultes(details = True, sort = True, only_identified = False)
#shot.plot( None, log_scale = True, sensitivity_correction = True)
#exit()
#shot.plot_all( None, log_scale = True, sensitivity_correction = True)
#exit()
#s = spectrometer.fittingFunction(array((1,0.1,0.5)),x)
#Shot( 1200, HR4000)
#print strptime('Thu Feb 04 16:27:51 CET 2010', "%a %b %d %H:%M:%S %Z %Y")
#shot = HR2000.loadData('./data_nova/4.02.10/', 'SpectraSuite_txt', 'HSA100610')
#figure()
#shot = HR2000.loadData('./data_nova/4.02.10/', 'SpectraSuite_txt', 'HSA10069')
#print strptime('Thu Feb 04 16:27:51 CET 2010', "%a %b %d %H:%M:%S %Z %Y")
#exit()
#shot = HR2000.loadData('./data_nova/10.10.2011/','SpectraSuite_xml','shot011')
#shot.shot_num = hash('HSA1011a0')
#print shot.shot_num
#spec = shot.spectra[0]
#spec.processSpectrum()
#shot.processShot()
#spec.plot()
#shot.saveData(noiseCorrected = True)
#shot.plot(None, log_scale = False, sensitivity_correction = True)
#shot.plot_all(None, log_scale = False, sensitivity_correction = False)
#shot.tableOfResultes(details = True, sort = True, only_identified = True)
#sh993a
#shot.processShot()
#shot.tableOfResultes(details = False, sort = True, only_identified = True)
#spectrometer.convertCountToPhotons(arange(200,1200),1,True)
#plot(x,s)
#show()
#exit()
#path = './data/'
#COMPASS = Tokamak('COMPASS')
#HR2000 = OceanOptics_Spec('HR+C1103',COMPASS)
#HR2000.integTime = 1.8#[ms]
#print COMPASS.actualShotNumber()
#HR2000.start()
#try:
#dir_list = os.listdir(path)
#for i in dir_list:
#os.remove(path+i)
#except OSError:
#os.mkdir(path)
#print 'can not remove old files'
#while True:
#print 'waiting on software trigger'
#COMPASS.waitOnSoftwareTrigger() #TODO otestovat
#print 'software trigger accepted'
#print 'waiting on hardware trigger'
##COMPASS.waitOnHardwareTrigger(path+'spectra.txt')
#print 'hardware trigger accepted'
#shot = HR2000.loadData(path,'Java_Data_file','spectra')
#shot.saveData(filename = None, noiseCorrected = True)
##shot.plot_all('png', log_scale = True, sensitivity_correction = True)
##shot.plot('png', log_scale = False, sensitivity_correction = True, w_interval = [200, 900])
#shot.tableOfResultes(details = False, sort = True, only_identified = False)
##shot.plotLinesTimeEvolution('png',only_identified = True)
#COMPASS.storeData(path)
#exit()
path = './data/'
GOLEM = Tokamak('GOLEM')
HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
HR2000.integTime = 1.62#[ms]
HR2000.start()
while True:
try:
dir_list = os.listdir(path)
for i in dir_list:
os.remove(path+i)
except OSError:
os.mkdir(path)
print 'can not remove old files'
print 'waiting on software trigger'
GOLEM.waitOnSoftwareTrigger() #TODO otestovat
print 'software trigger accepted'
print 'waiting on hardware trigger'
GOLEM.waitOnHardwareTrigger('ready')
print 'hardware trigger accepted'
shot = HR2000.loadData(path,'Java_Data_file','spectra')
shot.saveData(filename = None, noiseCorrected = True)
shot.plot_all('png', log_scale = True, sensitivity_correction = True)
shot.plot('png', log_scale = False, sensitivity_correction = True, w_interval = [200, 900])
shot.tableOfResultes(details = False, sort = True, only_identified = True)
shot.tableOfResultes(details = True,sufix = 'full', sort = True, only_identified = False)
shot.plotLinesTimeEvolution('png',only_identified = True)
GOLEM.storeData(path)
print 'data stored '+asctime(localtime())
print 'success !!!!!'
exit()
#exit()
#path = './data_nova/'
#GOLEM = Tokamak('GOLEM')
##HR4000 = OceanOptics_Spec('HR4C1947',GOLEM)
##HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
##HR2000 = OceanOptics_Spec('HR+C0732',GOLEM)
##HR2000 = OceanOptics_Spec('HR+C1834',GOLEM)
##HR2000 = OceanOptics_Spec('HR+C1103',GOLEM)
#working_dirs = os.listdir(path)
#database = list()
#for directory in working_dirs:
#shots = os.listdir(path+directory)
#for data_file in shots:
#if data_file[-8:] == 'SpectraSuite_xml':
#data_file = data_file[:-9]
#name = path+directory+'/','SpectraSuite_xml',data_file
#print name
#try:
#shot = HR2000.loadData(path+directory+'/','SpectraSuite_xml',data_file)
#shot.tableOfResultes(details = True, sort = True, only_identified = True)
#shot.plot_all(None, log_scale = False,sensitivity_correction = False)
#shot.saveData(filename = None, noiseCorrected = True)
#except Exception as inst:
#print 'problem', inst
path = './data_nova/'
COMPASS = Tokamak('COMPASS')
HR2000 = OceanOptics_Spec('HR+C0732',COMPASS)
working_dirs = os.listdir(path)
database = list()
integtime = list()
def ReadFile(directory,data_file):
name = path+directory+'/','txt',data_file
#print name
#if data_file[-3:] == 'txt':
#data_file = data_file[:-4]
#print name
#return
if data_file[-8:] == 'ProcSpec':
data_file = data_file[:-9]
#name = path+directory+'/','ProcSpec',data_file
#print name
try:
#print locale.getdefaultlocale()
shot = HR2000.loadData(path+directory+'/','SpectraSuite_xml',data_file)
#shot.tableOfResultes(details = True, sort = True, only_identified = True)
#if shot.readoutNoiseRMS != 0:
integtime.append(shot.integ_time)
if shot.integ_time > 6:
return
shot.plot_all(None, log_scale = True,sensitivity_correction = True)
#shot.spectra[0].plot()
#clf()
#show()
except Exception as inst:
print 'problem', inst
for directory in working_dirs:
shots = os.listdir(path+directory)
#if directory!= '2012':
#continue
for data_file in shots:
if os.path.isdir(path+directory+'/'+data_file):
shots_2 = os.listdir(path+directory+'/'+data_file)
for data_file_2 in shots_2:
#print path+directory+'/'+data_file+'/'+data_file_2
ReadFile(directory+'/'+data_file,data_file_2)
else:
ReadFile(directory,data_file)
#integtime = array(integtime)
#hist(integtime,50, range = (0, 20))
#show()
if __name__ == "__main__":
main()
# vrací to oto blbě výslednou tlouštku
#TODO otestovat zabranení vícenásobného spousteni
# dívat tu javu v nekonečné smyčce?
[Return]