#!/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/>.
  
"""



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))
	    
	
	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)):
			print 'can not copy icon'
		if os.path.isfile(data_folder+'status'):    
		    os.system(('scp -r %s %s:%s' % (data_folder+'status',self.login+'@'+self.uploadServerIp,self.database+'status'))%(self.shotNumber))os.system(('scp -r %s %s:%s' % (data_folder+'status',self.login+'@'+self.uploadServerIp,self.database+'status'))%(self.shotNumber))
		    print 'status ready'
		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' 
	#self.shotNumber = 7341
	#return    
	[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, 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 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)
	
	#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:
	    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()
	    #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(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))
	    #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 ./data/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 +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')
	    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 =  0#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]
	    
	    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)
	    
	    
	    
	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()
	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)
	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)
	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)
	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):
	"""
	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,  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#1.67#[ms]

    #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')
	clf()
	plot(shot.darkNoise)
	savefig(shot.DataFile+'_darkNoise.png')
	
	
	#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(5)
    
      
 
 
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 = './2900/'   
    
    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) 


    def ReadFile(directory):	
	
		
	# V. weinzettell format 

	try:
	    shot = HR2000.loadData(directory+'/','VW_java','')
	    shot.plot_all(None, log_scale = True,sensitivity_correction = False)
	    shot.plot(None, log_scale = False,sensitivity_correction = False)

	    shot.tableOfResultes(details = False, sort = True, only_identified = False)

    
	  
	    shot.plotLinesTimeEvolution(None,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:    
		#if './data_wladaserv/Halpha/2012_03_06/65' != path+s_dir+'/'+d_dir+'/'+data_file:
		    #continue

		print path+s_dir+'/'+d_dir+'/'+data_file
		ReadFile(path+s_dir+'/'+d_dir+'/'+data_file)
   

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):
    #plot a narrow interval around some line, compare different shots

     

    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():
     # interference filters accessible on GOMLEM tokamak
     
     
    GOLEM = Tokamak('GOLEM')
    HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
    
    shot_num =  8059
    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))    

    filters = loadtxt('golem_filters.txt')
    sensitivity_photodiod = loadtxt('photodiod_sensitivity.txt')
    sens_interp = interpolate.interp1d(sensitivity_photodiod[:,0],
    sensitivity_photodiod[:,1]/100,bounds_error=False, fill_value=0.0, kind = 'cubic')
    #plot(sens_interp(shot.wavelength), )
    #show()

    
    
    integ_spectrum = zeros(HR2000.resolution)
    for spectrum in shot.spectra:  
	if  spectrum.plasma():
	    integ_spectrum += spectrum.getData(sensitivity_correction = True)
    plot(shot.wavelength, integ_spectrum*sens_interp(shot.wavelength),linewidth=2)
    xlabel('wavelength [nm]')
    ylabel('counts [a.u.]')

    x = shot.wavelength
    for i,filtr in enumerate(filters):	
	plot( x,exp(-(x-filtr[0])**2/(2*(filtr[1]/2.35)**2))*10000)
	
    show()
    
   

 

    # MAIN PROGRAM
    
def main(): 

    Example1()

    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)


    #Example8()


    
    
    
if __name__ == "__main__":
    main()