#!/usr/bin/env python 
# -*- coding: utf-8 -*-
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
from  shutil import copyfile
import locale

from matplotlib.ticker import NullFormatter,ScalarFormatter

set_printoptions(precision=4,linewidth=400)

rcParams['backend'] = 'Agg'




def MAD(x):
    return median(abs(x-median(x)))*1.4826
    
def MovingAverage(data, subset_size):
    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):
    
    x = linspace(x[0],x[-1],size(x)*10) #increase resolution
    xc =p[2]
    A = p[0]
    w = 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
#TODO vše se tu počítá bez přístrojového rozšíření, u tisku parametrů dát možnost correct instrumental broadening

def lineCharacteristics(p):
    pc = zeros(5)
    a = p[3]
    n = p[4]

    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:
    def __init__(self, name):
	self.name = name
	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 compasu není N, v golemovi není B, i v pořadí v jakém se nají hledat

	
    def saveConfig(self):
	config = ConfigParser.RawConfigParser()
	config.set('DEFAULT', 'Name', self.name)
	
	config.add_section('Basic parameters')
	config.set('Basic parameters', 'Plasma length', self.dataAcqusitionTime)
	config.set('Database', 'remote database folder', self.database)
	config.set('Database', 'trigger port', self.trigger_port)
	config.set('Database', 'database login', self.login)

	
	with open('tokamak_'+self.name +'.cfg', 'wb') as configfile:
	    config.write(configfile)
 #GOLEM
 
    def storeData(self, data_folder):  
	
	#TODO odstranit!!!
	#self.uploadServerIp = '192.168.2.125'
	#self.shotNumber = 6879
	
	
	if self.uploadServerIp == '0.0.0.0' or self.shotNumber == -1:  #TODO opravit
	    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:		
		os.system(('scp -r %s %s:%s' % (data_folder,self.login+'@'+self.uploadServerIp,self.database ))%(self.shotNumber))
		if os.path.isfile(data_folder+'icon.png'):
		    os.system(('scp -r %s %s:%s' % (data_folder+'icon.png',self.login+'@'+self.uploadServerIp,self.database+'icon.png'))%(self.shotNumber))
    
		break
	    except Exception as inst:
		print 'problem with database access:', inst 
		sleep(5)
		
		
    def actualShotNumber(self):  #TODO, vrátí číslo aktual shot, u golema nastavené fcí wait trigger
	return self.shotNumber
	
    def waitOnSoftwareTrigger(self):  #čeká dokud nedorazí trigger, u golema to nastaví proměnnou aktual shot  
    #echo pokus | nc 192.168.0.12 50007 -q0

	def netcat(hostname, port):
	    s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
	    s.bind((hostname, port))
	    s.listen(1)
	    conn, addr = s.accept()
	    data = conn.recv(1024)
	    conn.close()
	    return addr,data
      
	HOST = ''                 # Symbolic name meaning all available interfaces
	PORT =  self.trigger_port            # Arbitrary non-privileged port
	[addr,data] = netcat(HOST, PORT)
	self.actualShotNumber = int(data)   # je to integer
	self.uploadServerIp = addr[0]
	
	
	
    def waitOnHardwareTrigger(self, filename):
	
	#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)
	    
	
	    
	    

    
        #os.stat(filename).st_mtime   # watch file creation time


class Spectrometer():
    "Spectrometer properties and methods container class " 
    def __init__(self, SerialNumber,tokamak):
	#self.spectrometerPrecision = 0
	#self.wavelength
	self.SerialNumber = SerialNumber
	self.Tokamak = tokamak
	
    #def setAcqusitionParameters(self,exposure = None, NumberOfSpectra)
    #def startDataAcqusition(self)
    #def peakFittingFunction(self,):

    #def dispersion(self,wavelength) #nm/px
    #def estimateNoiseInCounts(self,wavelength, intensity):
    #def convertCountToPhotons(self,wavelength, intensity, make_plot = False):
    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
    
    
    
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.elements = list()
	    #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')
	    for name in database:
		self.elements.append([loadtxt('SpectralLines/'+name),name[0:-4]])
	    
	    #load everything from config file
	    self.loadConfig()
	    
	    
	
    #def start(self):
	##os.system('trap "echo Koncim, zabijim javac ...;  killall javac ; exit"  SIGTERM SIGINT EXIT')
	#try:
	    #os.remove('nohup.out')
	    #os.remove('highspeedacquisitionsample/HighSpeedAcquisitionSample.class')
	#except:
	    #print 'file can not be removed' 
	
	#java_library = '/opt/OmniDriver-/OOI_HOME/HighResTiming.jar:/opt/OmniDriver-/OOI_HOME/OmniDriver.jar'
	
	#os.system('javac -classpath '+java_library+' HighSpeedAcquisitionSample.java')
	
	#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 nice -20 java -Djava.library.path=. -classpath '+java_library+':. highspeedacquisitionsample.HighSpeedAcquisitionSample %s %i %i &'  %run_options)

	#os.system('tail -f nohup.out &')
	
    def start(self):
	#os.system('trap "echo Koncim, zabijim javac ...;  killall javac ; exit"  SIGTERM SIGINT EXIT')
	try:
	    os.remove('nohup.out')
	    os.remove('highspeedacquisitionsample/HighSpeedAcquisitionSample_simple.class')
	except:
	    print 'file can not be removed' 
	
	java_library = '/opt/OmniDriver-/OOI_HOME/HighResTiming.jar:/opt/OmniDriver-/OOI_HOME/OmniDriver.jar'
	
	os.system('javac -classpath '+java_library+' HighSpeedAcquisitionSample_simple.java')
	
	copyfile('HighSpeedAcquisitionSample_simple.class', 'highspeedacquisitionsample/HighSpeedAcquisitionSample_simple.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('nice -20 java -Djava.library.path=. -classpath '+java_library+':. highspeedacquisitionsample.HighSpeedAcquisitionSample_simple %s %i %i &'  %run_options)

	os.system('tail -f nohup.out &')
	
	
	
	
		
	
	
    def wavelengthCorrection(self, old_wavelength):

	new_wavelegth = copy(old_wavelength)
	for i,a in enumerate(self.calibrationPolynomeCorrection):
	    new_wavelegth+= a*old_wavelength**i
	
	return  new_wavelegth
	
	
	
	
    def loadConfig(self):
	
	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_lenght = config.getfloat('Efficiencies', 'fiber lenght') 
	
	tmpstr = config.get('Other properties', 'Calibration polynome correction')
	self.calibrationPolynomeCorrection = 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):
	    
	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)  
	#TODO  Nonlinearity coefficients
	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 lenght', self.fiber_lenght)
	
	
	config.add_section('Other properties')
	config.set('Other properties', 'Calibration polynome correction', self.calibrationPolynomeCorrection)
	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):
	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] = sqrt(p[1]**2+self.F_inst_width(p[2])**2)	
	p_full[3] = self.F_inst_skewness(p[2])
	p_full[4] = self.voitig_parameter 
	
	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):	
	p_full = self.calcParamsForFitting(p)
	return lineCharacteristics(p_full)
    
	
    def elementIdentification(self,wavelength,sigma):

	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.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):

	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_lenght)
	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
	#wavelength = linspace(190,1000,1000)
	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):
	
	print 'peak fitting'
	yn = intensity[ix0]
	ixr = ix0+1 

	while ixr < self.resolution 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
	


	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

	popt_corr = self.lineCharacteristics(popt)
	
	chisq=sum(info["fvec"]*info["fvec"])
	doF=ixr-ixl-3
	#print "\t\tchi^2/doF ", chisq/max(doF,1), " limit ", scipy.stats.chi2.isf(0.1,doF)/max(doF,1)
	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]))

	colour = ( 'g','r','c', 'm', 'y', 'k' )
	if   line == inf:
	    plot_fit = True
	#plot data
	if plot_fit:
	    step(wavelength[max(ixl-3,0):min(ixr+3,self.resolution)],intensity[max(ixl-3,0):min(ixr+3,self.resolution)],'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,self.resolution))
	    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,self.resolution-1)],10*(ixr-ixl))
	    plot(interv,self.fittingFunction(popt,interv), '--')
	    axhline(y=4, ls= '--',label = 'threshold')
	    
	    #print('wavelength '+str(popt_corr[2])+' +/- '+str(sqrt(pcov[2,2])))
	    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") 

	    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')
	    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()


	return popt_corr,sqrt(diag(pcov)),ixl,ixr,line_name, line

    	
	
    def identifyLines(self,wavelength, intensity,readoutNoiseRMS,  debug = True):
	# return table of lines, wavelngths, +other parametres
	print '\t identifing lines'

	
	
	spectraParametres = list()
	line_fit = list()
	length = len(intensity)

	derivation = zeros(length)	
	
	mean_peak_FWHM = 5 #[px]  #TODO je to nesystematické, dát do vlastností spektrometru
    
	lim = scipy.stats.chi2.isf(0.01/length,mean_peak_FWHM)
	
	#intensity-= median(intensity[self.black_pixels])   #  dávat to sem? 
	#readoutNoiseRMS = MAD(intensity)
	intensity/=readoutNoiseRMS	
	
	
	while True:   

	    cumSumSpectra = cumsum((intensity*(intensity>0))**2)*2 - arange(length)
	    
	    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')
		#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')
		#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],(line-popt[2] ), sig[2]))

	    savetxt('current_peak.txt',(wavelength[ixl:ixr],intensity[ixl:ixr]))
	    # remove the fitted peak (in ideal case the peak should be substracted
	    intensity[max(ixl+1,0):min(ixr, self.resolution-1)+1] = 0

	return line_fit
	
	
	
    def loadData(self,path,typ,name = 'spectra'):
	
	
	shot = Shot(self.Tokamak.actualShotNumber(),  self)
	
	if typ == 'simple':   # 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])
		if line[:22] == 'Board temperature [C]:':
		    shot.temperature = float(line[24:-1])
		if line[:17] == 'Spectrometer S/N:':  #TODO opravit od do
		    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] #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'):		
		shot.darkNoise = genfromtxt(path+'background.txt', skip_header=5)
		shot.darkNoise -= median(shot.darkNoise[self.black_pixels])
	    
		f = open(path+'background.txt', 'r')
		for i,line in enumerate(f):
		    if i > 20:
			break
		    if line[0:22] == 'Integration time [us]:':
			integ_time_background = int(line[24:-1])
		shot.darkNoise*=shot.integ_time/float(integ_time_background)
		f.close()

	    
	    #extract data about readout patterns
	    if  os.path.isfile(path+'readoutpatterns.txt'):
		shot.readOutPatterns =  genfromtxt(path+'readoutpatterns.txt', skip_header=1)
		shot.readOutPatterns-= (shot.darkNoise/shot.integ_time)*self.min_integ_time*1000
		#shot.readOutPatterns -= median(shot.readOutPatterns[self.black_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
		#print 'shot.n_spectra', shot.n_spectra
		#print myzip.namelist()
		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 = zeros(shot.n_spectra)
	    #readOutPatterns = zeros(shot.n_pixels)
	    
	    #TODO odstranit
	    
	    #def hotPixels(DarkSpectra):
		#mad= MAD(DarkSpectra)
		#hotPixels = DarkSpectra - median(DarkSpectra)> +mad*3
		
		#plot(DarkSpectra)		
		#plot(where(hotPixels)[0],DarkSpectra[(hotPixels)],'ro')
		#show()
		#return where(hotPixels)[0] 
	    
	    #for i in range(4):
		#readOutPatterns+= loadtxt(path+'dark'+str(i+1)+'.txt')[:,1]
	    #readOutPatterns/= 4.0
	  
	    #TODO udělat inicializace shot.readoutNoiseRMS
	    #shot.readoutNoiseRMS  = MAD(intensity[self.black_pixels])
	    shot.readoutNoiseRMS = MAD(intensity)
	    print 'mad',  MAD(intensity[self.black_pixels]), MAD(intensity)
	    print 'integ time ', shot.integ_time 
	    #shot.readOutPatterns = MovingAverage(readOutPatterns,100)
	    
	    #self.hotPixels = hotPixels(readOutPatterns-shot.readOutPatterns)   

	    #plot(intensity)
	    #plot((0,2000), (shot.readoutNoiseRMS*4,4*shot.readoutNoiseRMS))
	    #show()
	    
	    #shot.readOutPatterns -= median(shot.readOutPatterns)  
	   
	   
	   
	if typ == 'SpectraSuite_txt':		#text acqurid by Data Studio
	    #TODO tested only for 1 acqusitioned spectra per file!
	    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())
	    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 wavelangth 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 = ''
	
	
	#TODO dopsat sem všechny položky, která vytváří load data
	
	# pamatuje si pole objektů Spektrum
	
    # ve složce path musí být soubory ,spectra.txt,background.txt,readoutpatterns.txt
   
    
    
    def plasmaShot(self):  #TODO zjistit jestli tam bylo plazma    
	pass

    def getData(self,sensitivity_correction = False,noiseCorrected = True): 
	
	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):  # vrátí první sloupec vlnové délkym zbylé 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/1000.0)+'\n'
	header += 'Acquisition time [ms]:\t'+	'%2.2f' %(self.time_stamps[-1]/(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):
	pass   #TODO dodělat
	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):
	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)

	
	for i,spectrum in enumerate(self.spectra):
	    ax = fig.add_subplot(self.n_spectra ,1,i+1)
	    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(self.n_spectra ,1,self.n_spectra )  
	ax.xaxis.set_major_formatter( ScalarFormatter() )
	
	
	if file_type == None:
	    show()
	else:
	    savefig(self.DataFile+'_Graph_.'+file_type)
	close()

    def processShot(self, details = False):
	
	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):
	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-1,0,-1):
	    if self.spectra[i].plasma():
		plasma_end = i+2
		break 
	
	print IntensitiesTable
	IntensMax =  amax(IntensitiesTable, axis = 1)

	sort_ind = argsort(IntensMax, order=None)[::-1]	

	#zipped =  zip(NamesTable,IntensitiesTable,ErrorsTable)	


	    
	for i,index in enumerate(sort_ind):
	    if i > 6:
		break
	    errorbar(time_stamps[:plasma_end],IntensitiesTable[index][:plasma_end],yerr=ErrorsTable[index][:plasma_end],label=NamesTable[index])
	    
	legend(loc = 'best')
	axis([0,time_stamps[plasma_end-1],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 tableOfResultes(self, sort = False, details = False, only_identified = True):
	
	ColoumnNames,LinesTable,time_stamps = self.processShot(details) 
	
	f = open(self.DataFile+'_Resultes.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():
    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):
	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):
	
	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):

	intensity = copy(self.intensity)
	
	if noiseCorrected:
	    intensity -= self.Shot.readOutPatterns+self.Shot.darkNoise
	    intensity[int_(self.Shot.Spectrometer.hotPixels)] = median(intensity)
	    
	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)	
	self.lineFitList = self.Shot.Spectrometer.identifyLines(self.wavelength,corr_intensity,self.Shot.readoutNoiseRMS,  False ) # return table of lines, wavelngths, +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 *
    #GOLEM = Tokamak('GOLEM')
    #Spectrometer('gfdg',GOLEM)
    #HR4000 = OceanOptics_Spec('HR4C1947',GOLEM)
    #HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)

    
    ##HR2000 = OceanOptics_Spec('HR+C1834',GOLEM)

    #HR2000 = OceanOptics_Spec('HR+C1103',GOLEM)
    

    #shot = HR2000.loadData('./data_nova/20.08.2010/', 'SpectraSuite_xml', 'shot000') 
    
    #shot.plot_all(None)
    
    #shot.tableOfResultes(details = True, sort = True, only_identified = True)

    #shot = HR2000.loadData('./data_golem/','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') 
    #HR2000 = OceanOptics_Spec('HR+C0732',GOLEM)  
    #shot = HR2000.loadData('./data_nova/02.12.2011/','SpectraSuite_xml','shot005') 
    #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/'
    GOLEM = Tokamak('GOLEM')
    HR2000 = OceanOptics_Spec('HR+C1886',GOLEM)
    
    HR2000.integTime = 2.0#[ms]
    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'
	#GOLEM.waitOnSoftwareTrigger()			#TODO otestovat
	#print 'software trigger accepted'
	
	#print 'waiting on hardware trigger'	
	#GOLEM.waitOnHardwareTrigger(path+'spectra.txt')
	#print 'hardware trigger accepted'
	
	#shot = HR2000.loadData(path,'simple','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.plotLinesTimeEvolution('png',only_identified = True)
	#GOLEM.storeData(path) 
	
    #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/'   
    
    GOLEM = Tokamak('GOLEM')

    HR2000 = OceanOptics_Spec('HR+C0732',GOLEM)  
    
    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()
    	
    	
    	
    	
    	
    
    	
# opravit     def storeData(self, data_folder):  
 #otesovat software trigger
 
# jak zjistit koeficienty nelinearity?