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]

Navigation