source: ETALON/CLIO/control/read_plot_data.py @ 744

# -*- coding: utf-8 -*-
"""
Created on Mon Feb  6 15:47:25 2017

@author: delerue
"""


import numpy as np
import matplotlib.pyplot as plt
#from time import sleep
from shutil import copyfile,move
import os.path
from scipy.fftpack import fft, ifft
from scipy import signal

def read_plot_data(filename):
    #Definitions
    nchannels=32
        
    #electron_signal
    electrons_channel=26
    electrons_signal_scaling=0.02
    electrons_threshold=-1000
    electrons_threshold_peak=200.
    electrons_half_width=1

    data_start=1
    data_length=8000
    electron_pretrig=5   
#    electron_trig_stop=500
    electron_trig_stop=20
    electron_post_trig_integration_start=1
    electron_post_trig_integration_stop=500

    #position signal
    position_channel=28

    #channels mapping
    #channels_mapping=[17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]
    data_mapping=[ 0,   111,   83,   76,   48, 
                   55,  104,   69,   62,    0,
                   97,   90,    0,  118,    0,
                   0,   125,    0,    0,    0,
                   0,     0,    0,    0,    0,
                   0,     0,    0,    0,    0,
                   0,   0 ]


    fft_depth=150#lower value -- lower number of freq.

    #figure_size (in inches, that is 100*pixels with dpi=100)
    figure_size=(8, 6)
    figure_dpi=100
    
    #sort the angles and find corresponding channels
    angles_sorted=[]
    idx_angles_sorted=[]
    idx_data_mapping=np.zeros(len(data_mapping),dtype=np.int)
    
    for angle in np.sort(data_mapping):
        if angle>0:
            if len(angles_sorted)==0:
                angles_sorted=[angle]
                idx_angles_sorted=[data_mapping.index(angle)]
            else:
                angles_sorted.append(angle)
                idx_angles_sorted.append(data_mapping.index(angle))
            idx_data_mapping[data_mapping.index(angle)]=int(len(angles_sorted))
    #Check if the filename contains the path or not
    if (len(os.path.dirname(filename))==0):
        filename="./"+filename
    dayname=os.path.basename(filename).replace("-","").split("_")[2]
    target_directory=os.path.dirname(filename)+"/../"+dayname+"/"
    target_filename=target_directory+os.path.basename(filename)
    
    #read parameters    
    key_values={}
    with open(filename,'r') as fid:
        fileline=fid.readline()
        while (len(fileline)>0)and(fileline[0]=='#'):                
            entries=fileline.replace('#','').replace(' ','').replace('\n','').split('=')
            key_values[entries[0]]=entries[1]
            fileline=fid.readline()
    fid.close()
        
    for p in key_values: print p,'=', key_values[p]       
    #read data    
    alldata=np.loadtxt(filename, dtype='int32', comments='#', delimiter='  ', usecols=range(0,nchannels), converters={ _:lambda s:  int(s , 16)  for _ in range(0,nchannels) })
    
    #converting the data in the correct format
    alldata=alldata%(2**16)-2**15;

        #Scaling
    #data_scaling_factor=(float(key_values['amplitude'])*5.)/(2.**15)
    data_scaling_factor=1;
    time_scaling_factor=(float(key_values['acquisition_rate'])/36.)*1.e-3

        #remove constant components from electron signal
    FFT_spectrum=fft(alldata[1:data_length,electrons_channel]);
    Amplitude=abs(FFT_spectrum);
    Amplitude[0]=0;#remove constant component
    alldata[1:data_length,electrons_channel]=np.real(ifft(Amplitude*np.exp(1j*np.angle(FFT_spectrum))));


        #look for the electrons
    elec=alldata[1:data_length,electrons_channel]
    pos_elec=[];
    Elec_sig_val_arr=[];
    step=round(data_length/(2*(data_length*time_scaling_factor/40.)))#40 ms repetition rate
    for c in range(1,int(round(data_length/step)+1)):
        llim=int((c-1)*step);
        hlim=int(step*c);
        if hlim>data_length:
            hlim=int(data_length)
        if (max(elec[llim:hlim])>electrons_threshold_peak):
            pos_elec.append(llim+np.argmax(elec[llim:hlim]))
            Elec_sig_val_arr.append(min(elec[llim:hlim]));
    #print Elec_sig_val_arr

    electrons_found=False

        # fft filtering
    for c in range(0,11):
                FFT_spectrum=fft(alldata[1:data_length,idx_angles_sorted[c]]);
                Phase=np.angle(FFT_spectrum);
                Amplitude=abs(FFT_spectrum);
                Amplitude[(1+fft_depth):(data_length-fft_depth)]=0;
                Amplitude[0]=0;#remove constant component
                Signal=ifft(Amplitude*np.exp(1j*Phase));
                alldata[1:data_length,idx_angles_sorted[c]]=np.real(Signal);

        #Find the amplitude of the signal
    Sig = [0 for x in range(len(pos_elec))] 
    Sig_mean=[0 for x in range(12)]
    Sig_std=[0 for x in range(12)]
    Elec_sig_val=0;
    Elec_sig_val_err=0;
    if (len(pos_elec)>0.):
        electrons_found=True
        if (pos_elec[0]-50)>0:
            electron_data_start=pos_elec[0]-50
        else:
            electron_data_start=0
        if (pos_elec[0]+50)<data_length:
            electron_data_stop=pos_elec[0]+50
        else:
            electron_data_stop=(data_length-1)        
        Elec_sig_val=np.mean(Elec_sig_val_arr);
        Elec_sig_val_err=np.std(Elec_sig_val_arr);
        print 'electrons found:', Elec_sig_val,'+/-', Elec_sig_val_err
                # Signal amplitude
        for c in range(0,12):
            for k in range(0,len(pos_elec)):
                Sig[k]=alldata[(pos_elec[k]+14),idx_angles_sorted[c]]-alldata[(pos_elec[k]-14),idx_angles_sorted[c]]
            Sig_mean[c]=-np.mean(Sig)#mean of signal in one file
            Sig_std[c]=np.std(Sig)#std in one file

        move(filename,target_filename)
            
    #print Sig_mean, Sig_std

 
    #prepare the figures    
    fig=plt.figure(1,figsize=figure_size,dpi=figure_dpi)
    if electrons_found:
        figzoom=plt.figure(2,figsize=figure_size,dpi=figure_dpi)
        figsignal=plt.figure(3,figsize=figure_size,dpi=figure_dpi)

        #set channels to plot
    plotch=data_mapping;
    plotch[electrons_channel]=1;
    plotch[16]=0;
   # plotch[electrons_channel+1]=1001;

    for idata in range(0,nchannels-1):
        
        this_data_scaling_factor=data_scaling_factor
        txtline='r'
        linestyle='-';
        if data_mapping[idata] >0:
            txtline='r'
            if data_mapping[idata] >1000:
                txtline='b'
                
            if idata==electrons_channel:
                        txtline='k'
                        this_data_scaling_factor=this_data_scaling_factor*electrons_signal_scaling
            plt.figure(1)
            plt.plot(np.arange(data_start,data_length)*time_scaling_factor,alldata[range(data_start,data_length),idata]*this_data_scaling_factor, txtline+linestyle)
            if (electrons_found):
                plt.figure(2)                       
                plt.plot(np.arange(electron_data_start,electron_data_stop)*time_scaling_factor, (alldata[range(electron_data_start,electron_data_stop),idata]-alldata[(pos_elec[0]-14),idata])*this_data_scaling_factor, txtline+linestyle)

            
    plt.figure(1)
    plt.title(key_values['date']+ ' ---  Motor value:  '+key_values['motor_value']+ ' ---  Current:  '+key_values['current_value'])
    plt.ylabel('Signal amplitude [a.u.]')
    plt.xlabel('Time [ms]')
    plt.grid(True)
    imagename=target_filename.replace('.txt','.png')
    print(imagename)
    fig.savefig(imagename)
    plt.close(fig)
    imagename_no_path=os.path.basename(imagename)


    fid=open(target_directory+'list.html','a');
    fid.write("<A HREF=")
    fid.write(imagename_no_path)
    fid.write("><IMG SRC=")
    fid.write(imagename_no_path)
    fid.write(" width=600>")
    fid.write("</A><BR/>")
    fid.write(imagename_no_path)
    fid.write("<BR/><BR/>\n")
    fid.close()
    copyfile(imagename,os.path.dirname(imagename)+'/last.png')

    if (electrons_found):
        plt.figure(2)
        plt.title(key_values['date']+ ' ---  Motor value:  '+key_values['motor_value']+ ' ---  Current:  '+key_values['current_value'])
        plt.ylabel('Signal amplitude [a.u] - pedestal removed')
        plt.xlabel('Time [ms]')
        plt.grid(True)        
        imagenamezoom=target_filename.replace('.txt','_zoom.png')
        figzoom.savefig(imagenamezoom)
        plt.close(figzoom)
        copyfile(imagenamezoom,os.path.dirname(imagename)+'/last_zoom.png')
    
        imagenamezoom_no_path=os.path.basename(imagenamezoom)
        fid=open('zoom_list.html','a');
        fid.write("<A HREF=")
        fid.write(imagenamezoom_no_path)
        fid.write("><IMG SRC=")
        fid.write(imagenamezoom_no_path)
        fid.write(" width=600>")
        fid.write("</A><BR/>")
        fid.write(imagenamezoom_no_path)
        fid.write("<BR/><BR/>\n")
        fid.close()
        
        plt.figure(3)
        plt.plot(range(48,132,7), Sig_mean,'r^-')
        plt.title(key_values['date']+ ' ---  Motor value:  '+key_values['motor_value']+ ' ---  Current:  '+key_values['current_value'])
        plt.ylabel('Signal amplitude [a.u.]')
        plt.xlabel('Detector angle [deg.]')
        plt.grid(True)        
        imagenamesignal=target_filename.replace('.txt','_signal.png')
        figsignal.savefig(imagenamesignal)
        plt.close(figsignal)
        copyfile(imagenamesignal,os.path.dirname(imagename)+'/last_signal.png')

        imagenamesignal_no_path=os.path.basename(imagenamesignal)
        fid=open('signal_list.html','a');
        fid.write("<A HREF=")
        fid.write(imagenamesignal_no_path)
        fid.write("><IMG SRC=")
        fid.write(imagenamesignal_no_path)
        fid.write(" width=600>")
        fid.write("</A><BR/>")
        fid.write(imagenamesignal_no_path)
        fid.write("<BR/><BR/>\n")
        fid.close()
        
    filenamesig=target_filename.replace('.txt','.sig')
    fid=open(filenamesig,'w');
    for key in key_values:
            fid.write("#" +key+" = "+key_values[key]+" \n");
    fid.write("#electronsig = "+ str(Elec_sig_val)+" \n");  
    fid.write("#electronstd = "+ str(Elec_sig_val_err)+" \n");      
    for idata in range(0,len(Sig_mean)):
        fid.write(str(angles_sorted[idata])+"\t"+str(Sig_mean[idata])+"\t"+str(Sig_std[idata])+"\n")
    fid.write("\n")
    fid.close()

    filenamedone=filename.replace('.txt','.done')
    fid=open(filenamedone,'a');
    fid.write("done")
    fid.close()
    
    return;