# -*- coding: utf-8 -*-
"""
Created on Tue Mar 22 15:23:25 2016

@author: Martian
"""
import matplotlib.pyplot as plt
from numpy import genfromtxt,exp,sqrt,diag
import numpy as np
import scipy as sp
from scipy.optimize import curve_fit


def ProcessData(filename):
    with open(filename,"r") as data:
        for line in data:
            if line.startswith("&END"):
                break
        else:
            print "Oooops,failed to read data"
        stored_data=genfromtxt(data,skip_header=1)

        raw_data=np.array(stored_data,dtype=float)
        ch1_5keV=raw_data[raw_data[:,1]==1,0]
        ch2_5keV=raw_data[raw_data[:,1]==2,0]
        ch1_10keV=raw_data[raw_data[:,3]==1,2]
        ch2_10keV=raw_data[raw_data[:,3]==2,2]
        ch1_15keV=raw_data[raw_data[:,5]==1,4]
        ch2_15keV=raw_data[raw_data[:,5]==2,4]
        ch1_20keV=raw_data[raw_data[:,7]==1,6]
        ch2_20keV=raw_data[raw_data[:,7]==2,6]
        ch1_25keV=raw_data[raw_data[:,9]==1,8]
        ch2_25keV=raw_data[raw_data[:,9]==2,8]

        return (ch1_5keV,ch2_5keV,ch1_10keV,ch2_10keV,ch1_15keV,ch2_15keV,ch1_20keV,ch2_20keV,ch1_25keV,ch2_25keV)


def Histograms(ch1_10keV,ch2_10keV):

    t=np.linspace(0,10,500)
    t2=(t[:-1]+t[1:])/2.0
    plt.figure()
    plt.hist(ch1_10keV,bins=t)
    plt.title("py14mz Channel Left 10 keV")
    plt.xlabel("Time(us)")
    plt.ylabel("Counts(s-1)")

    plt.figure()
    plt.hist(ch2_10keV,bins=t)
    plt.title("py14mz Channel Right 10 keV")
    plt.xlabel("Time(us)")
    plt.ylabel("Channel Counts(s-1)")
    plt.show()

    return (t2)

def Probabilities(ch1_5keV,ch2_5keV,ch1_10keV,ch2_10keV,ch1_15keV,ch2_15keV,ch1_20keV,ch2_20keV,ch1_25keV,ch2_25keV):

    t=np.linspace(0,10,500)

    counts_left_5keV=np.array((plt.hist(ch1_5keV,bins=t))[0],dtype=float)
    counts_right_5keV=np.array((plt.hist(ch2_5keV,bins=t))[0],dtype=float)
    counts_left_10keV=np.array((plt.hist(ch1_10keV,bins=t))[0],dtype=float)
    counts_right_10keV=np.array((plt.hist(ch2_10keV,bins=t))[0],dtype=float)
    counts_left_15keV=np.array((plt.hist(ch1_15keV,bins=t))[0],dtype=float)
    counts_right_15keV=np.array((plt.hist(ch2_15keV,bins=t))[0],dtype=float)
    counts_left_20keV=np.array((plt.hist(ch1_20keV,bins=t))[0],dtype=float)
    counts_right_20keV=np.array((plt.hist(ch2_20keV,bins=t))[0],dtype=float)
    counts_left_25keV=np.array((plt.hist(ch1_25keV,bins=t))[0],dtype=float)
    counts_right_25keV=np.array((plt.hist(ch2_25keV,bins=t))[0],dtype=float)

    P_L_5keV=(counts_left_5keV)/(counts_right_10keV.sum()+counts_left_10keV.sum())
    P_R_5keV=(counts_right_5keV)/(counts_right_10keV.sum()+counts_left_10keV.sum())
    P_L_10keV=(counts_left_10keV)/(counts_right_10keV.sum()+counts_left_10keV.sum())
    P_R_10keV=(counts_right_10keV)/(counts_right_10keV.sum()+counts_left_10keV.sum())
    P_L_15keV=(counts_left_15keV)/(counts_right_10keV.sum()+counts_left_10keV.sum())
    P_R_15keV=(counts_right_15keV)/(counts_right_10keV.sum()+counts_left_10keV.sum())
    P_L_20keV=(counts_left_20keV)/(counts_right_10keV.sum()+counts_left_10keV.sum())
    P_R_20keV=(counts_right_20keV)/(counts_right_10keV.sum()+counts_left_10keV.sum())
    P_L_25keV=(counts_left_25keV)/(counts_right_10keV.sum()+counts_left_10keV.sum())
    P_R_25keV=(counts_right_25keV)/(counts_right_10keV.sum()+counts_left_10keV.sum())

    return (counts_left_5keV,counts_right_5keV,counts_left_10keV,counts_right_10keV,counts_left_15keV,counts_right_15keV,counts_left_20keV,counts_right_20keV,counts_left_25keV,counts_right_25keV,P_L_5keV,P_R_5keV,P_L_10keV,P_R_10keV,P_L_15keV,P_R_15keV,P_L_20keV,P_R_20keV,P_L_25keV,P_R_25keV,)


def Assymetry_data(P_L_5keV, P_R_5keV,P_L_10keV,P_R_10keV,P_L_15keV,P_R_15keV,P_L_20keV,P_R_20keV,P_L_25keV,P_R_25keV):

    A_5keV=(P_L_5keV-P_R_5keV)/(P_L_5keV+P_R_5keV)
    A_10keV=(P_L_10keV-P_R_10keV)/(P_L_10keV+P_R_10keV)
    A_15keV=(P_L_15keV-P_R_15keV)/(P_L_15keV+P_R_15keV)
    A_20keV=(P_L_20keV-P_R_20keV)/(P_L_20keV+P_R_20keV)
    A_25keV=(P_L_25keV-P_R_25keV)/(P_L_25keV+P_R_25keV)

    return (A_5keV,A_10keV,A_15keV,A_20keV,A_25keV)


def GuessB (A_5keV,A_10keV,A_15keV,A_20keV,A_25keV):
        B_10keV=np.argmax(A_10keV)

        return(B_10keV)


def A (t2, beta,B):

    gamma=851.616
    return (-1.0/3.0)*((sp.sin(gamma*B*t2 - beta)-sp.sin(gamma*B*t2 + beta))/(2.0*beta))

def AssymetryErr (counts_left_5keV,counts_right_5keV,counts_left_10keV,counts_right_10keV,counts_left_15keV,counts_right_15keV,counts_left_20keV,counts_right_20keV,counts_left_25keV,counts_right_25keV):

    delta_right_5keV=sqrt(counts_right_5keV)
    delta_left_5keV=sqrt(counts_left_5keV)
    delta_right_10keV=sqrt(counts_right_10keV)
    delta_left_10keV=sqrt(counts_left_10keV)
    delta_right_15keV=sqrt(counts_right_15keV)
    delta_left_15keV=sqrt(counts_left_15keV)
    delta_right_20keV=sqrt(counts_right_20keV)
    delta_left_20keV=sqrt(counts_left_20keV)
    delta_right_25keV=sqrt(counts_right_25keV)
    delta_left_25keV=sqrt(counts_left_25keV)

    err_5keV = 2*sqrt(1/(counts_left_5keV+counts_right_5keV)**4*(counts_left_5keV**2*delta_right_5keV**2+counts_right_5keV**2*delta_left_5keV**2))
    err_10keV = 2*sqrt(1/(counts_left_10keV+counts_right_10keV)**4*(counts_left_10keV**2*delta_right_10keV**2+counts_right_10keV**2*delta_left_10keV**2))
    err_15keV = 2*sqrt(1/(counts_left_15keV+counts_right_15keV)**4*(counts_left_15keV**2*delta_right_15keV**2+counts_right_15keV**2*delta_left_15keV**2))
    err_20keV = 2*sqrt(1/(counts_left_20keV+counts_right_20keV)**4*(counts_left_20keV**2*delta_right_20keV**2+counts_right_20keV**2*delta_left_20keV**2))
    err_25keV = 2*sqrt(1/(counts_left_25keV+counts_right_25keV)**4*(counts_left_25keV**2*delta_right_25keV**2+counts_right_25keV**2*delta_left_25keV**2))

    return (err_5keV,err_10keV,err_15keV,err_20keV,err_25keV)

def Assymetry_Plot(t2,A_5keV,A_10keV,A_15keV,A_20keV,A_25keV,err_5keV,err_10keV,err_15keV,err_20keV,err_25keV):

    plt.figure()
    plt.plot(t2,A_10keV)
    plt.title("py14mz Asymmetry Function 10 keV")
    plt.xlabel("Time(us)")
    plt.ylabel("Channel Counts(s-1)")
    plt.show

    fitParams_5keV, fitCovariance_5keV = curve_fit(A,t2,A_5keV,po=(1.116,0.008),absolute_sigma=True,sigma=err_5keV)
    fitParams_10keV, fitCovariance_10keV = curve_fit(A,t2,A_10keV,po=(1.116,0.008),absolute_sigma=True,sigma=err_10keV)
    fitParams_15keV, fitCovariance_15keV = curve_fit(A,t2,A_15keV,po=(1.116,0.008),absolute_sigma=True,sigma=err_15keV)
    fitParams_20keV, fitCovariance_20keV = curve_fit(A,t2,A_20keV,po=(1.116,0.008),absolute_sigma=True,sigma=err_20keV)
    fitParams_25keV, fitCovariance_25keV = curve_fit(A,t2,A_25keV,po=(1.116,0.008),absolute_sigma=True,sigma=err_25keV)

    plt.plot(t2, A(t2,fitParams_10keV[0],fitParams_10keV[1]), 'r-', label='fit')
    plt.errorbar(t2,A_10keV,yerr=err_10keV,fmt='',color='g',alpha=0.4)

    return (fitParams_5keV,fitParams_10keV,fitParams_15keV,fitParams_20keV,fitParams_25keV)

ProcessData("pr3.dat")
#Histograms(ch1_10keV,ch2_10keV)
Probabilities(ch1_5keV,ch2_5keV,ch1_10keV,ch2_10keV,ch1_15keV,ch2_15keV,ch1_20keV,ch2_20keV,ch1_25keV,ch2_25keV)
Assymetry_data(P_L_5keV,P_R_5keV,P_L_10keV,P_R_10keV,P_L_15keV,P_R_15keV,P_L_20keV,P_R_20keV,P_L_25keV,P_R_25keV)
AssymetryErr(counts_left_5keV,counts_right_5keV,counts_left_10keV,counts_right_10keV,counts_left_15keV,counts_right_15keV,counts_left_20keV,counts_right_20keV,counts_left_25keV,counts_right_25keV)
Assymetry_Plot(A_10keV)