Source code :: SpectrometerControl
[Return]
[Download]#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Ocean Optics spectrometer data analyzer
Program for controlling of the Ocean optics spectrometers and extractin the most important information from spectra
## oospecanalyz is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 2 of the License, or
## (at your option) any later version.
## oospecanalyz is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
## You should have received a copy of the GNU General Public License
## along with oospecanalyz. If not, see <http://www.gnu.org/licenses/>.
f
"""
from numpy import *
import matplotlib
matplotlib.use('agg')
from matplotlib.pyplot import *
import os
from time import mktime, strptime,sleep,asctime, localtime,struct_time, 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
import urllib2
from matplotlib.ticker import NullFormatter,ScalarFormatter
set_printoptions(precision=4,linewidth=400)
rcParams['backend'] = 'Agg'
debugging = False
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_interp = interpolate.interp1d(arange(len(x)),x)
x0 = x
m = 10
x = linspace(x[0],x[-1]+(x[-1]-x[-2]),size(x)*m)-mean(diff(x0))/2#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(y)/m,m) )
y = sum(y,1)/m #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])
#http://en.wikipedia.org/w/index.php?title=Skew_normal_distribution
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') # list of the element - originaly comes from NIST database
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
"""
if self.shotNumber is None:
print 'missing shot number'
return
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))
for i in range(10):# neměla by tou ta smyšby být
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)):
print 'can not copy icon'
if os.path.isfile(data_folder+'status'):
if os.system(('scp -r %s %s:%s' % (data_folder+'status' ,self.login+'@'+self.uploadServerIp,self.database+'status' ))%(self.shotNumber)):
print 'can not copy'
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
[addr,data] = netcat(HOST, PORT)
if data == '':
self.shotNumber = None
print 'wrong shot number reclieved'
else:
self.shotNumber = int(data) # je to integer
self.uploadServerIp = addr[0]
def waitOnHardwareTrigger(self, filename, t ):
"""
Detect that spectrometer has reclieved hardware trigger
"""
#detect the hardware trigger only by change of the file with the data
t0 = time()
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 and time()-t0 < t :
if not os.path.isfile(filename):
change_time = 0
else:
change_time = os.stat(filename).st_mtime
sleep(0.1)
if creat_time != change_time:
print 'hw trigger ', str(time()-t0 ), 's'
return True
else:
print 'timeout'
return False
#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
if not isnan(sigma):
presicion = min(sqrt(self.wavelength_presicion**2+sigma**2),dispersion*3)
else:
presicion = min(sqrt(self.wavelength_presicion**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,integ_time, 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.
"""
if self.absolute_calibration == None and self.relative_calibration == None:
f1 = interpolate.interp1d(self.sensor_efficiency[:,0],self.sensor_efficiency[:,1], bounds_error=True, fill_value=0.3, kind = 'quadratic')
f2 = interpolate.interp1d(self.gratting_efficiency[:,0],self.gratting_efficiency[:,1], bounds_error=True, fill_value=0.1, kind = 'quadratic')
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, kind = 'quadratic')
intensity/= (f1(wavelength)*f2(wavelength)*f3(wavelength))*self.photonsPerCount
if make_plot:
plot(wavelength,self.photonsPerCount/(f1(wavelength)*f2(wavelength*f3(wavelength))))
else:
if self.absolute_calibration:
f = interpolate.interp1d( self.absolute_calibration[:,0], self.absolute_calibration[:,1], kind = 'slinear')
intensity/= f(wavelength)
else:
f = interpolate.interp1d( self.relative_calibration[:,0], self.relative_calibration[:,1], kind = 'slinear')
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, kind = 'quadratic')
intensity/= f(wavelength)*f3(wavelength)
#intensity/= integ_time/1000.0
if make_plot:
plot(wavelength,f(wavelength))
if make_plot:
xlabel('wavelength [nm]')
ylabel('unit/count')
show()
return intensity
def peakFitting(self,ix0,wavelength,intensity, plot_fit = False):
"""
Identify the peak around the middle in index ix0. 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)
x0 = wavelength[ix0]
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)
#print popt_corr
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
#print line
#if abs(line - 657.805 )<0.01:
#plot_fit = True
if popt_corr[1] > 5:
popt_corr[0] = nan
#plot_fit = True
#plot graph for easier debugging
if plot_fit:
clf()
interv = arange(max(ixl-3,0),min(ixr+3,n))
w_interp = interpolate.interp1d(interv,wavelength[interv])
crop_w = wavelength[interv]
step(crop_w,intensity[interv],'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
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
crop_w_100 = w_interp(linspace(amin(interv), amax(interv), len(interv)*100))
crop_w_1 = w_interp(linspace(amin(interv), amax(interv), len(interv)*1))
plot(crop_w_100,self.fittingFunction(popt,crop_w_100), 'b--')
plot(wavelength[ixl:ixr],self.fittingFunction(popt,wavelength[ixl:ixr]), '.')
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")
#CM = sum(wavelength[ixl:ixr]*self.fittingFunction(popt,wavelength[ixl:ixr]))/sum(self.fittingFunction(popt,wavelength[ixl:ixr]))
#axvline(x = CM, ls = '-.',c = 'y', label = "Center of mass 2")
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, A: %4.1f'%(popt_corr[2],sqrt(pcov[2,2]),chisq/doF,popt_corr[0] ))
show()
#savefig(str(self.Tokamak.shotNumber)+'.png')
#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')
xlim(amin(wavelength),amax(wavelength))
show()
popt,sig,ixl,ixr,line_name,line = self.peakFitting(ix0,wavelength,intensity,debug)
#renormalize
popt[0] = self.convertCountToPhotons(x0,readoutNoiseRMS*popt[0],1 )
sig[0] = self.convertCountToPhotons(x0,readoutNoiseRMS*sig[0] ,1)
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))
#print line_fit
#exit()
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('./data/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 > ./data/nohup_java.out &' %run_options)
#os.system('nohup chrt -f 99 java -Djava.library.path=. -classpath '+java_library+':. highspeedacquisitionsample.HighSpeedAcquisitionSample %s %i %i &' %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)
try:
tmpstr = config.get('Efficiencies', 'Absolute calibration') # counts/(W*m-2*nm)
self.absolute_calibration = parseStringOfField(tmpstr)
except ConfigParser.NoOptionError as (err):
self.absolute_calibration = None
try:
tmpstr = config.get('Efficiencies', 'Overall effectivity')
self.relative_calibration = parseStringOfField(tmpstr)
except ConfigParser.NoOptionError as (err):
self.relative_calibration = None
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 +self.min_integ_time#shift ->0
clf()
plot(diff(shot.time_stamps),'.', label = 'exposure times from timestamps')
axhline(y=self.integTime, ls= '--',label = ' exposure time ')
ylim(0,None)
legend(loc = 'best')
savefig(shot.DataFile+'_diff_timestamps.png',bbox_inches='tight')
clf()
print 'read out time', mean(diff(shot.time_stamps)), ' ; ', diff(shot.time_stamps)
# I think, that timestams returned by HighSpeedAcqusion are not the real timestamps of the spectra
shot.time_stamps = self.Tokamak.trigger_shift+self.min_integ_time*2+(arange(shot.n_spectra)+0.5)*shot.integ_time
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:])/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
if line[0:22] == 'Integration time [us]:':
integ_time_readoutpatterns = int(line[24:-1])/1000
if line[0:16] == 'RMS^2 of signal:':
shot.readoutNoiseRMS = sqrt(float(line[17:]))
print 'RMS', shot.readoutNoiseRMS
f.close()
shot.readOutPatterns = genfromtxt(path+'readoutpatterns.txt', skip_header=imax)
if abs(mean(shot.readOutPatterns[:1023]) -mean(shot.readOutPatterns[1023:])) > 5: # SPECTROMETER FAILURE
shot.readOutPatterns[1023:] -= shot.readOutPatterns[1024]-shot.readOutPatterns[1023]
shot.readoutNoiseRMS = 10
shot.readOutPatterns-= (shot.darkNoise/shot.integ_time)*integ_time_readoutpatterns
#plot(shot.readOutPatterns)
#show()
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')
line_header = 0
for i,line in enumerate(f):
line_header +=1
ind = line.find('Serial Number')
if ind != -1:
SN = line[ind+15:-1]
ind = line.find('Date and time (GMT)')
if ind != -1:
#print (line[ind+21:-1])
#print strptime("Thu Apr 5 21:46:30 2012", "%a %b %d %H:%M:%S %Y") #BUG!!!!
#print strptime((line[ind+21:-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:
shot.readoutNoiseRMS = float(line[ind+20:-1])
if line.find('***************') != -1:
break
if SN != self.SerialNumber:
raise Exception('wrong serial number, '+SN+' != '+self.SerialNumber)
f.close()
if shot.readoutNoiseRMS > 10:#probably failura in the calculation
shot.readoutNoiseRMS = 10
f = open(shot.DataFile, 'r')
data = genfromtxt(f,skip_header=line_header)
f.close()
shot.wavelength = data[:,0]
#print shot.n_spectra
#print shape(data[:,i+1])
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 == 'VW_java': #acqurid by V.W. java program
shot.temperature = nan
shot.DataFile = path
if not os.path.isfile(shot.DataFile+'wavelength_axis.txt'):
raise Exception('Data file '+shot.DataFile+'wavelength_axis.txt'+' does not exists.')
shot.wavelength = genfromtxt(shot.DataFile+'wavelength_axis.txt')
if max(abs(nanmin(shot.wavelength)-self.wavelength_range[0]), abs(nanmax(shot.wavelength)-self.wavelength_range[1]))>10:
raise Exception('bad wavelength range %d - %d != %d - %d' %(nanmin(shot.wavelength),nanmax(shot.wavelength),self.wavelength_range[0],self.wavelength_range[1]))
if not os.path.isfile(shot.DataFile+'settings.txt'):
raise Exception('Data file '+shot.DataFile+'settings.txt'+' does not exists.')
f = open(shot.DataFile+'settings.txt', 'r')
for i,line in enumerate(f):
#Number of measured spectra:
if i == 1:
shot.n_spectra = int(line)
#Requested integration time [microseconds]:
if i == 3:
shot.integ_time = float(line)
f.close()
shot.time = os.stat(shot.DataFile+'settings.txt').st_mtime
shot.n_pixels = 2048
if not os.path.isfile(shot.DataFile+'timing.txt'):
raise Exception('Data file '+shot.DataFile+'timing.txt'+' does not exists.')
times = genfromtxt(shot.DataFile+'timing.txt',skip_header=2)
shot.time_stamps = times[:,0]-times[0,0] +self.Tokamak.trigger_shift #[ms]
if os.path.isfile(shot.DataFile+'temperature.txt'):
shot.temperature = genfromtxt(shot.DataFile+'temperature.txt')
for i in arange(shot.n_spectra):
if not os.path.isfile(shot.DataFile+'spectrum_'+str(i+1)+'.txt'):
raise Exception('Data file '+shot.DataFile+'spectrum_'+str(i+1)+'.txt'+' does not exists.')
intensity = loadtxt(shot.DataFile+'spectrum_'+str(i+1)+'.txt')
intensity -= median(intensity)
shot.spectra.append(Spectrum(shot,shot.wavelength,intensity, shot.time_stamps[i]))
shot.readoutNoiseRMS = MAD(intensity)
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.')
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.')
replaceComma = lambda s:float(s.replace(',','.'))
if os.path.isfile(path+'readoutpatterns.txt'):
readOutPatterns = genfromtxt(path+'readoutpatterns.txt',skip_header=17,skip_footer = 1,converters = {0:replaceComma, 1:replaceComma })
shot.readOutPatterns = readOutPatterns[:,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')
n_aver = 1
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[:16] == 'Spectra Averaged':
if line[-11:-1] == '('+self.SerialNumber+')':
n_aver = int(line[17:-11])
else:
n_aver = int(line[17:-1])
if line[:5] == '>>>>>':
break
if shot.n_spectra == 0: # missing information in file
shot.n_spectra = 1
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))
if n_aver == 1:
shot.readoutNoiseRMS = MAD(intensity)
else:
shot.readoutNoiseRMS = 10/sqrt(n_aver)# std(intensity[self.black_pixels])
intensity -= median(intensity[self.black_pixels])
#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(nanmin(shot.wavelength)-self.wavelength_range[0]), abs(nanmax(shot.wavelength)-self.wavelength_range[1]))>10:
raise Exception('bad wavelength range %d - %d != %d - %d' %(nanmin(shot.wavelength),nanmax(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 '
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, linearityCorrected = False)
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('units [-]')
else:
ylabel('counts [-]')
if file_type == None:
show()
else:
savefig(self.DataFile+'_Graph_%d.'%(i)+file_type,bbox_inches='tight')
close()
clf()
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,bbox_inches='tight')
close()
#clf()
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:
# every detect except the firts, is false
if isinstance(row[2+time_slice],float) and row[2+time_slice] != 0:
break
if isinstance(row[2+time_slice],(list,tuple)) and row[2+time_slice][0] != 0:
break
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+1)+'-th spectrum'
for line in Table: # lines are sorted from most intense
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 = 7, logscale = False):
"""
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
plasma_end = min(plasma_end,self.n_spectra-1)
plasma_start = 0
for i in arange(self.n_spectra):
if self.spectra[i].plasma():
plasma_start = i-1
break
plasma_start = max(0,plasma_start)
t_start = time_stamps[plasma_start]
t_end = time_stamps[plasma_end]
IntensMax = nanmax(IntensitiesTable, axis = 1)
sort_ind = argsort(IntensMax, order=None)[::-1]
for i,index in enumerate(sort_ind):
if i >= n_lines:
break
points= [IntensitiesTable[index] != 0]
if not logscale:
points[0][plasma_end] = True
points[0][plasma_start] = True
errorbar(time_stamps[points],IntensitiesTable[index][points],yerr=ErrorsTable[index][points],label=NamesTable[index])
leg = legend(loc='best', fancybox=True)
leg.get_frame().set_alpha(0.5)
axis( [t_start,t_end,0, nanmax(IntensitiesTable) *1.1])
xlabel('t [ms]')
ylabel('Intensity [a.u.]')
title('The most intense lines')
if logscale:
yscale('log', nonposy='clip')
ylim([nanmin(IntensitiesTable)*0.8, nanmax(IntensitiesTable)*1.1])
if file_type == None:
show()
else:
savefig(self.DataFile+'_LinesEvolution_.'+file_type,bbox_inches='tight')
close()
clf()
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 = [nanmin(self.wavelength), nanmax(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,self.Shot.integ_time)
else:
noise *= ones(self.Shot.Spectrometer.resolution)
min_threshold = nanmin(threshold*noise)
noise *= w_interval_bool
max_threshold = nanmax(threshold*noise)
plot(self.wavelength, intensity,'b', label= '%2.1f ms' %self.time_stamp, linewidth=0.3)
plot(self.wavelength,threshold*noise, 'g', linewidth=0.6)
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,nanmax(intensity))])
else:
axis([w_interval[0],w_interval[1],nanmin(-noise) , max(max_threshold,nanmax(intensity))])
axis([w_interval[0],w_interval[1],0, max(max_threshold,nanmax(intensity))])
leg = legend(loc='center right', fancybox=True)
leg.get_frame().set_alpha(0.7)
def getData(self,sensitivity_correction = False,noiseCorrected = True, linearityCorrected = 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)
if linearityCorrected:
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,self.Shot.integ_time)
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, debugging ) # return table of lines, wavelength, +other parametres
return self.lineFitList
def Example0():
# the most primitive example of accesing data from GOLEM - store them to the same folder as this program, name of file with data is spectra.txt_Data_.txt
GOLEM = Tokamak('GOLEM')
HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
shot_num = 8420
GOLEM.shotNumber = shot_num
import urllib2
url = 'http://golem.fjfi.cvut.cz/operation/shots/'+str(shot_num)+'/diagnostics/Radiation/1111Spectrometer.ON/data/spectra.txt_Data_.txt'
u = urllib2.urlopen(url)
localFile = open('data_'+str(shot_num)+'.txt', 'w')
localFile.write(u.read())
localFile.close()
shot = HR2000.loadData('','Python_Data_file','data_'+str(shot_num))
print shot.wavelength
for spectrum in shot.spectra:
print spectrum.intensity
shot.plot(None, log_scale = False, sensitivity_correction = True)
shot.plot_all(None, log_scale = False, sensitivity_correction = True)
shot.plotLinesTimeEvolution(None,only_identified = True)
def Example1():
#run acqusition on Golem
path = './data/'
GOLEM = Tokamak('GOLEM')
HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
HR2000.integTime = 1.67#[ms] it doesn't work waster on this computer
#os.system('killall java')
#sleep(0.5)
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'
os.system('killall java')
sleep(0.5)
#software trigger - cca 20s before discharge
print 'waiting on software trigger'
GOLEM.waitOnSoftwareTrigger()
print 'software trigger accepted'
#prepare spectrometer (starting and measuring of the bacground take about 10s)
HR2000.start()
if not GOLEM.waitOnHardwareTrigger('prepared', 3): # very problematics step!!!!!
print os.system('killall java')
print 'java killed'
sleep(0.5)
HR2000.start()
#hardware trigger - watche if the drivers hase stored the date
print 'waiting on hardware trigger'
if not GOLEM.waitOnHardwareTrigger('ready',60):
print 'hardware trigger timeout'
continue
print 'hardware trigger accepted'
shot = HR2000.loadData(path,'Java_Data_file','spectra')
shot.saveData(filename = None, noiseCorrected = True)
#ploting data
shot.plot_all('svgz', log_scale = True, sensitivity_correction = True)
shot.plot('svgz', log_scale = False, sensitivity_correction = True, w_interval = [200, 900])
shot.plot_all('png', log_scale = True, sensitivity_correction = True)
shot.plot('png', log_scale = False, sensitivity_correction = True, w_interval = [200, 900])
plot(shot.readOutPatterns)
savefig(shot.DataFile+'_readoutnoise.png',bbox_inches='tight')
clf()
plot(shot.darkNoise)
savefig(shot.DataFile+'_darkNoise.png',bbox_inches='tight')
#processing spectra
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)
shot.plotLinesTimeEvolution('svgz',only_identified = True)
#create status file
File = open('./data/status', 'w')
File.write('0') #status READY
File.close()
GOLEM.storeData(path)
print 'data stored '+asctime(localtime())
sleep(10)
def Example2():
# load file "spectra.txt_Data_" (file from internet, where was extracted all information about spectra) and horizontal lines at the position of the spectral lines at database
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('None',only_identified = True)
colour = ( 'g','r','c', 'm', 'y', 'k' )
#plot positions of all ions at to the spectra together
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()
#plot positions of all ions separately
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()
def Example3():
# explore primitive database of COMPASS spectra, data must be stored in folder './data_wladaserv/'
path = '../../data_wladaserv/'
COMPASS = Tokamak('COMPASS')
#HR2000 = OceanOptics_Spec('HR+C0732',COMPASS) #VIS
HR2000 = OceanOptics_Spec('HR+C1834',COMPASS) #H alpha
#HR2000 = OceanOptics_Spec('HR+C1103',COMPASS) #UV
spectrometer_dirs = os.listdir(path)
print spectrometer_dirs
def ReadFile(directory):
# V. weinzettell format
try:
shot = HR2000.loadData(directory+'/','VW_java','')
shot.plot_all('png', log_scale = True,sensitivity_correction = False)
#shot.plot_all('png', log_scale = True,sensitivity_correction = False)
shot.plot(None, log_scale = False,sensitivity_correction = True)
#shot.tableOfResultes(details = False, sort = True, only_identified = False)
shot.plotLinesTimeEvolution('pdf',only_identified = True, logscale = True)
shot.plotLinesTimeEvolution('png',only_identified = True, logscale = True)
except Exception as inst:
print 'problem', inst
for s_dir in spectrometer_dirs:
date_dirs= os.listdir(path+s_dir)
for d_dir in date_dirs:
shots_dir= os.listdir(path+s_dir+'/'+d_dir)
for data_file in shots_dir:
dir_path = path+s_dir+'/'+d_dir+'/'+data_file
print dir_path
#if os.path.isfile(dir_path+'/_Graph_.png'):
#continue
ReadFile(dir_path)
def Example4():
#explore on COMPASS "database", process data in ProcSpec files
path = './data_nova/'
GOLEM = Tokamak('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:] == '.ProcSpec':
data_file = data_file[:-9]
name = path+directory+'/','.ProcSpec',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 = False)
except Exception as inst:
print 'problem', inst
def Example5():
#run the acqusition on the COMPASS
path = './data/'
COMPASS = Tokamak('COMPASS')
HR2000 = OceanOptics_Spec('HR+C1103',COMPASS)
HR2000.integTime = 15#[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() #condensators charge???
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)
shot.tableOfResultes(details = False, sort = True, only_identified = False)
shot.plotLinesTimeEvolution('png',only_identified = True)
COMPASS.storeData(path)
def Example6():
#calibration data
path = './calibration_8-3-2012/'
GOLEM = Tokamak('GOLEM')
HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
shot = HR2000.loadData(path,'SpectraSuite_txt','Hg2')
shot.plot(None, log_scale = False)
w, s = shot.getData(sensitivity_correction = True)
print sum(s)
for i, element in enumerate(GOLEM.elements):
shot.spectra[0].plot( log_scale = False, sensitivity_correction = True)
for j,line in enumerate(element[0]):
axvline(x = line , ls = '--', c = 'r')
title(element[1])
ylabel('counts [a.u.]')
show()
def Example7():
#plot spectra and positions of the lines, it can be useful for identification of the lines
shots = [7818,]+range(8059,8075)
He_gas = arange(8059,8070)
H_gas = arange(8070,8075)
GOLEM = Tokamak('GOLEM')
HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
for shot_num in shots:
gas = ''
if any(He_gas == shot_num) :
gas = 'He'
if any(H_gas == shot_num) :
gas = 'H'
print 'shot number '+str(shot_num)
GOLEM.shotNumber = shot_num
#download data
url = 'http://golem.fjfi.cvut.cz/operation/shots/'+str(shot_num)+'/diagnostics/Radiation/1111Spectrometer.ON/data/spectra.txt_Data_.txt'
try:
u = urllib2.urlopen(url)
except:
print 'shot number '+str(shot_num)+ ' does not exist'
continue
localFile = open('data_'+str(shot_num)+'.txt', 'w')
localFile.write(u.read())
localFile.close()
#read data
shot = HR2000.loadData('','Python_Data_file','data_'+str(shot_num))
integ_spectrum = zeros(HR2000.resolution)
for spectrum in shot.spectra:
if spectrum.plasma():
integ_spectrum += spectrum.intensity
#plot positions of the ions lins
colour = ( 'g','r','c', 'm', 'y', 'k' )
for i, element in enumerate(GOLEM.elements):
for j,line in enumerate(element[0]):
axvline(x = line , c = colour[(i%len(colour))])
title('ion: '+element[1]+'; shot #'+str(shot_num)+'; fill gas: ' + gas)
plot(shot.wavelength,integ_spectrum)
show()
def Example8(shots, names):
#Assisted line identifications
GOLEM = Tokamak('GOLEM')
HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
#download data
for shot_num in shots:
if os.path.isfile('data_'+str(shot_num)+'.txt'):
continue
print 'shot number '+str(shot_num)
GOLEM.shotNumber = shot_num
url = 'http://golem.fjfi.cvut.cz/operation/shots/'+str(shot_num)+'/diagnostics/Radiation/1111Spectrometer.ON/data/spectra.txt_Data_.txt'
try:
u = urllib2.urlopen(url)
except:
print 'shot number '+str(shot_num)+ ' does not contain spectrometer data'
continue
localFile = open('data_'+str(shot_num)+'.txt', 'w')
localFile.write(u.read())
localFile.close()
#plot spectra
while(True):
print 'set the wavelength interval for the spectral lines comparison'
l = raw_input('from [nm] (200 default): ')
if l == '':
l = 200
l = double(l)
r = raw_input('to [nm] (1200 default): ')
if r == '':
r = 1200
r = double(r)
shots_list = list()
for shot_num in shots:
shots_list.append(HR2000.loadData('','Python_Data_file','data_'+str(shot_num)))
integ_spectrum = list()
for shot in shots_list:
integ_spectrum.append(zeros(HR2000.resolution))
for spectrum in shot.spectra:
if spectrum.plasma():
integ_spectrum[-1] += spectrum.getData()
wavelength = shots_list[0].wavelength
interval = where((wavelength> l) * (wavelength < r))
colour = ( 'g','r','k', 'm', 'b', 'c' )
style = ('-', '--', '-.', ':')
shift = 0.4
for i,shot_spectrum in enumerate(zip(integ_spectrum,names)):
plot(wavelength[interval],shot_spectrum[0][interval],style[i%4],linewidth=2, label = shot_spectrum[1])
#plot lines from database
for i, element in enumerate(GOLEM.elements):
axvline(x = 0 , label = element[1], c = colour[(i/len(style))], ls = style[i%len(style)])
for j,line in enumerate(element[0]):
if line> r or line < l:
continue
axvline(x = line+shift , c = colour[(i/len(style))], ls = style[i%len(style)])
text(line+shift , 0, element[1])
leg = legend(loc='center right', fancybox=True)
leg.get_frame().set_alpha(0.7)
xlim(l, r)
show()
def Example9():
#massive data procesing
GOLEM = Tokamak('GOLEM')
HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
#i = %
#download data
for shot_num in arange(9000,9200):
if not os.path.isfile('data_'+str(shot_num)+'.txt'):
#continue
print 'shot number '+str(shot_num)
GOLEM.shotNumber = shot_num
url = 'http://golem.fjfi.cvut.cz/operation/shots/'+str(shot_num)+'/diagnostics/Radiation/1111Spectrometer.ON/data/spectra.txt_Data_.txt'
try:
u = urllib2.urlopen(url)
except:
#print 'shot number '+str(shot_num)+ ' does not contain spectrometer data'
continue
localFile = open('data_'+str(shot_num)+'.txt', 'w')
localFile.write(u.read())
localFile.close()
shot = HR2000.loadData('','Python_Data_file','data_'+str(shot_num))
#print shot.spectra[1])
i+=1
#if i%10!= 0:
for spectrum in shot.spectra:
if spectrum.plasma():
ha = spectrum.getData()[1008:1013]
plot( ha/mean(ha))
ylim(0,None)
show()
#plot(shot.spectra[1])
#show()
#exit()
##plot spectra
#while(True):
##print 'set the wavelength interval for the spectral lines comparison'
##l = raw_input('from [nm] (200 default): ')
##if l == '':
##l = 200
##l = double(l)
##r = raw_input('to [nm] (1200 default): ')
##if r == '':
##r = 1200
##r = double(r)
##shots_list = list()
##for shot_num in shots:
##shots_list.append(
#shot = HR2000.loadData('','Python_Data_file','data_'+str(shot_num))#)
#shot.
##integ_spectrum = list()
##for shot in shots_list:
##integ_spectrum.append(zeros(HR2000.resolution))
##for spectrum in shot.spectra:
##if spectrum.plasma():
##integ_spectrum[-1] += spectrum.getData()
##wavelength = shots_list[0].wavelength
##interval = where((wavelength> l) * (wavelength < r))
##colour = ( 'g','r','k', 'm', 'b', 'c' )
##style = ('-', '--', '-.', ':')
##shift = 0.4
##for i,shot_spectrum in enumerate(zip(integ_spectrum,names)):
##plot(wavelength[interval],shot_spectrum[0][interval],style[i%4],linewidth=2, label = shot_spectrum[1])
###plot lines from database
##for i, element in enumerate(GOLEM.elements):
##axvline(x = 0 , label = element[1], c = colour[(i/len(style))], ls = style[i%len(style)])
##for j,line in enumerate(element[0]):
##if line> r or line < l:
##continue
##axvline(x = line+shift , c = colour[(i/len(style))], ls = style[i%len(style)])
##text(line+shift , 0, element[1])
##leg = legend(loc='center right', fancybox=True)
##leg.get_frame().set_alpha(0.7)
##xlim(l, r)
##show()
# MAIN PROGRAM
def Example10():
#Assisted line identifications
#exit()
GOLEM = Tokamak('GOLEM')
HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
#i = %
N = list()
center = list()
data = list()
data_full = list()
pressure = list()
temperature = list()
dates = list()
#download data
sesion_lines = zeros((0,28))
shot_names = list()
PlasmaLengths = list()
path = './net_data/'
wrong_shot = [7342,7732,7737,7312,7400,7580,7582,7583]
for shot_num in arange(5000,9700):
print shot_num
if shot_num in wrong_shot:
print 'wrong'
continue
#print os.path.isfile(path+'data_'+str(shot_num)+'.txt')
#print os.path.isfile(path+'p_'+str(shot_num)+'.txt')
#exit()
if not os.path.isfile(path+'data_'+str(shot_num)+'.txt') or not os.path.isfile(path+'p_'+str(shot_num)+'.txt') or not os.path.isfile(path+'s_'+str(shot_num)+'.txt') or not os.path.isfile(path+'e_'+str(shot_num)+'.txt') :
#continue
if os.path.isfile(path+'no_data_'+str(shot_num)):
print 'neni'
continue
#print 'loadim'
#shot_num = 9685
url = 'http://golem.fjfi.cvut.cz/operation/shots/'+str(shot_num)+'/diagnostics/Radiation/1111Spectrometer.ON/data/spectra.txt_Data_.txt'
url_p = 'http://golem.fjfi.cvut.cz/operation/shots/'+str(shot_num)+'/Initial_PfeifferMerkaVakua'
url_s = 'http://golem.fjfi.cvut.cz/operation/shots/'+str(shot_num)+'/basicdiagn/PlasmaStart'
url_e = 'http://golem.fjfi.cvut.cz/operation/shots/'+str(shot_num)+'/basicdiagn/PlasmaEnd'
try:
u = urllib2.urlopen(url)
u_p = urllib2.urlopen(url_p)
u_s = urllib2.urlopen(url_s)
u_e = urllib2.urlopen(url_e)
except:
localFile = open(path+'no_data_'+str(shot_num), 'w')
localFile.close()
print 'neni'
continue
localFile = open(path+'data_'+str(shot_num)+'.txt', 'w')
localFile.write(u.read())
localFile.close()
localFile = open(path+'p_'+str(shot_num)+'.txt', 'w')
localFile.write(u_p.read())
localFile.close()
localFile = open(path+'s_'+str(shot_num)+'.txt', 'w')
localFile.write(u_s.read())
localFile.close()
localFile = open(path+'e_'+str(shot_num)+'.txt', 'w')
localFile.write(u_e.read())
localFile.close()
print 'download'
#exit()
shot = HR2000.loadData('','Python_Data_file',path+'data_'+str(shot_num))
if not shot.plasmaShot():
continue
#print shot.spectra[1])
#i+=1
#if i%10!= 0:
print 'shot number '+str(shot_num)
GOLEM.shotNumber = shot_num
wavelength, intensities = shot.getData()
#print shape(intensities)
for spectrum in shot.spectra:
if spectrum.plasma():
#if any(spectrum.getData()[-1] >5000):
#plot(spectrum.getData())
#show()
data.append(single(copy(spectrum.getData())))
temperature.append(shot.temperature)
pressure.append(loadtxt(path+'p_'+str(shot_num)+'.txt'))
#print shot.time
dates.append(shot.time)
#plot(spectrum.getData())
#show()
#shot.plot(None, log_scale = False, sensitivity_correction = True, w_interval = [200, 900])
#show()
#intens = sum(intensities, 1)
#data_full.append(intens)
#print spectra
#CII = [284.0,427.1,658.5,678,712]
#CIII = [230.1]
#HeI = [389.4,587.5,668,706.7,729.6]
#HeII = [542.7]
#H = [486.1,656.4]
#NII = [399.6,489.9,500.2,568.2,617,648.4]
#OII = [375.6,435,441.7,539.2,664.2,672.4,777.6,845]
#ions = [CII, CIII, H, HeI, HeII, NII,OII]
#shot_lines = zeros(0)
shot_names.append(shot_num)
#pressure.append(loadtxt(path+'p_'+str(shot_num)+'.txt'))
length = float(loadtxt(path+'e_'+str(shot_num)+'.txt'))-float(loadtxt(path+'s_'+str(shot_num)+'.txt'))
PlasmaLengths.append(length)
#print length
#exit()
#TODO to heII kecá!!
#TODO dělit délkou plazmatu
#for i,ion in enumerate(ions):
#lines = zeros(len(ion))
#for j,line in enumerate(ion):
#interval = where(abs(wavelength- line)<0.7)
#lines[j] = sum(intens[interval])
##shot_lines = hstack((shot_lines,lines))
#ion = H
#lines = zeros(len(ion))
#for j,line in enumerate(ion):
#interval = where(abs(wavelength- line)<0.7)
#lines[j] = sum(array(intens[interval]) )
#N.append(sum(lines)/length)
#interval = where(abs(wavelength- 468.1)<0.7)[0]
#center.append( sum(interval* intens[interval]/sum(intens[interval]) ))
##axvline(x = line-0.5)
##axvline(x = line+0.5)
##plot(wavelength[interval], intens[interval])
##title(str(i)+ ' '+str(line) )
##show()
##for spectrum in shot.spectra:
##if spectrum.plasma():
##wavelength =
##ha = spectrum.getData()[1008:1013]
##plot( ha/mean(ha))
##ylim(0,None)
##show()
##print shape(sesion_lines), shape(shot_lines)
##sesion_lines = vstack((sesion_lines, shot_lines))
#center = array(center)
#shot_names = array(shot_names)
#pressure = array( pressure)
#temperature = array(temperature)
##print shape(array( pressure))
#N = (array(N))
##plot(shot_names, temperature, '.')
##show()
##TODO shift/teplota
#data = array(data)
#print shape(data)
###print data
#data_full = array(data_full)
##exit()
#ix0 = argmin(abs(778-wavelength))
#ix0 = 1011
###plot()
#inter = arange(198,206)
#inter = arange(198,300)
##for i in arange(size(data_full,0)):
##plot( data_full[i,:])
##show()
#Hbeta = zeros((size(data_full,0)))
#Herror = zeros((size(data_full,0)))
#for i in arange(size(data_full,0)):
##Tokamak.shotNumber
#GOLEM.shotNumber = shot_names[i]
#ix0 = 1011
##ix0 = argmax(data_full[i,(ix0-5):(ix0+5)])
##plot(wavelength[(ix0-50):(ix0+50)], data_full[i,(ix0-50):(ix0+50)])
##show()
#ix0 = argmax(data_full[i,(ix0-3):(ix0+3)])+ix0-3
#popt,sig,ixl,ixr,line_name,line = HR2000.peakFitting(ix0,wavelength,data_full[i,:], plot_fit = False)
#Hbeta[i] = popt[2]
#Herror[i] = sig[2]
##print popt[2], sig[2]
##Hbeta[i,:] = [popt[2], sig[2]]
#print shape(arange(size(data_full,0))), shape(Hbeta)
##plot(arange(size(data,0)), Hbeta[:,0])
##print type(Hbeta[:,0])
##print type(Hbeta[:,1])
##print shape(Hbeta)
##print
##plot(arange(size(data,0)), Hbeta[:,1],'.')
#plot(shot_names,-Hbeta,'.')
##errorbar( Hbeta, Herror,'.')
#show()
#plot(temperature,Hbeta)
##errorbar( Hbeta, Herror,'.')
#show()
#plot(shot_names,temperature)
##errorbar( Hbeta, Herror,'.')
#show()
#plot(data.T)
#show()
save('data',data)
save('wavelength',wavelength)
save('temperature', temperature)
save('pressure', pressure)
save('shot_names', shot_names)
save('dates', dates)
#exit()
##print (N), (shot_names)
##plot(shot_names, N, '.')
#plot(pressure, N, '.')
#xlabel('$p_{init}$ [mPa]')
#ylabel('intensity H I')
#ylim(0,None)
#show()
#lengths = array(lengths)
#plot(pressure, lengths, '.')
#show()
#plot(shot_names, lengths, '.')
#show()
#plot(shot_names, center)
#show()
#print shape(shot_names)
#plot(shot_names, pressure, '.')
#ylabel('$p_{init}$ [mPa]')
#xlabel('shot number')
#ylim(0,None)
#show()
#sesion_lines = sesion_lines.T
##print shape(sesion_linessesion_lines)
#savetxt('fyztyd_lines',sesion_lines, fmt='%.2e' )
#print sesion_lines
#print array(shot_names)
#savetxt('fyztyd_lines_names',shot_names )
def main():
Example1()
#Example10()
#exit()
#shots = raw_input("Enter discharge numbers seperated by a comma (default 7818,8059,8074): ")
#if shots == '':
#shots = [7818, 8059, 8074]
#names = ('H, baked chamber','He, dirty chamber','H, dirty chamber')
#else:
#names = [ (i) for i in shots.replace(' ', '').split(',') ]
#shots = int_(names) #strip the string of spaces and split numbers based on comma and pack into a list of ints
###names = str_(shots)
###shots = [ int(i) for i in shots.replace(' ', '').split(',') ]
##print names, size(names), shots, size(shots)
##exit()
#Example8(shots, names)
if __name__ == "__main__":
main()[Return]
