figure()
sttl = 'Background - R1017'
cts = msc349bkg
ttl = 'IP=1000ms, 20 scan summation, H2counts={:.0f}, H2/H20={:.2f}%'
ttl = ttl.format(cts[1], (cts[1]/cts[14])*100)
plot_MS_bar(massBkg, cts, sttl, ttl)

#==========Generating cross talk plots from data====================#

#specific case for 2015-349 containing data for both H2 and water
#the script did not terminate D003 correctly
#raw D003 contains H2 OSNB CS-off, H2O bkgc, H2O bkgo, H2O OSNB CS-On
ds349 = get_datasets(sorted(msd349.keys()))
ds349['D003']['srange'] = (101, 140)
ds349['D000bkgc'] = ds349.pop('D000')
_ = ds349.pop('D004')
_ = ds349.pop('D005')

ds349h2o = {'D000bkgc': {'srange': (151, 170)},
            'D000bkgo': {'srange': (171, 190)},
            'D003': {'srange': (196, 235)},
            'D004': {'srange': (236, 275)},
            'D005': {'srange': (276, 315)}}

#ds349 uses the same OSNB background dataset as ds349h2o
ds349['D000bkgo']['sum'] = ds349h2o['D000bkgo'].copy()

#include 'H2O' in the ds349h2o sttl
for k in ds349h2o.keys():
    ds349h2o[k]['sttl'] = ds349h2o[k]['sttl']+' H2O'
    print ds349h2o[k]['sttl']

#the background scans in ds349 are masses 1-50 -> truncate to 1-23
dst = ds349
for k in dst.keys():
    if len(dst[k]['masses']) > 19:
        dst[k]['masses'] = masses
        dst[k]['sum'] = dst[k]['sum'][:19]

#the background scans in ds349h2o are masses 1-50 -> truncate to 1-23
dst = ds349h2o
for k in dst.keys():
    if len(dst[k]['masses']) > 19:
        dst[k]['masses'] = masses
        dst[k]['sum'] = dst[k]['sum'][:19]


#specific case for ds356
ds356['D004']['srange'] = (121, 145)
ds356['D001bkgc'] = ds356.pop('D001')
ds356['D002bkgo'] = ds356.pop('D002')

#specific case for ds015
ds015['D003']['sttl'] = 'OSNB - OS Off - R1023'
ds015['D001bkgc'] = ds015.pop('D001')
ds015['D002bkgo'] = ds015.pop('D002')

#specific case for ds018
ds018['D001bkgc'] = ds018.pop('D001')
ds018['D002bkgo'] = ds018.pop('D002')

#specific case for ds027
ds027['D001bkgc'] = ds027.pop('D001')
ds027['D002bkgo'] = ds027.pop('D002')

#specific case for ds028
_ = ds028.pop('D000')
ds028['D001bkgc'] = ds028.pop('D001')
ds028['D002bkgo'] = ds028.pop('D002')

#special case for ds046
#some of the scans in D000 occur between D002 and D003.
#Therefore need to sort on integer value of scan number
msc046.sort(key=lambda ms: int(ms[2]))
_ = ds046.pop('D000')
ds046['D001bkgc'] = ds046.pop('D001')
ds046['D002bkgo'] = ds046.pop('D002')

#specific case for ds047
#Observed evolution of outgassing products as filament warms up
#R1031-D003was moved to 2016-047 directory from 2016-048
#plot evolution of mass species

#specific case for ds048
#must be sorted similar to ds046
msc048.sort(key=lambda ms: int(ms[2]))
_ = ds048.pop('D000')
ds048['D004bkgc'] = ds048.pop('D004')
ds048['D005bkgo'] = ds048.pop('D005')


#specific case for 2015-352
dsets = get_datasets(msf)
msc352 = get_INMS_scan_data(msd352, id0='mass', sfields=sfld, mfields=mfld)
ds352 = {k:{'srange':(int(dsets[k][0][1:]),int(dsets[k][-1][1:]))}
			for k in dsets.keys()[1:4]}
ds352['D000bkgo'] = {'srange':(57,96)}
ds352['D000bkgc'] = {'srange':(17,56)}
_ = ds352.pop('D000')

plot_crosstalk(ds352, msc352)
dkeys = sorted(ds352.keys())
dsum = [ds352[k]['sum'] for k in dkeys]
dsum[2] = dsum[2]-dsum[1]
dsum[3] = dsum[3]-dsum[0]
dsum[4] = dsum[4]-dsum[1]

for i in range(2,5):
    sttl = ds352[dkeys[i]]['sttl']+' Bkg Subtract'
    cts = dsum[i]
    ttl = ds352[dkeys[i]]['ttl']
    plot_MS_bar(masses, cts, sttl, ttl)

cstoff = dsum[4]/dsum[3]*100
sttl = 'Crosstalk Ratio, 8.112kms, CS Off - R1020'
cts = cstoff
ttl = 'IP=1000ms, 40 scan summation, KE=0.675eV, H2 ratio={:.3f}%'.format(cts[1])
plot_MS_bar(masses, cts, sttl, ttl, yl='Percent')

cston = dsum[2]/dsum[3]*100
sttl = 'Crosstalk Ratio, 8.112kms, CS On - R1020'
cts = cston
ttl = 'IP=1000ms, 40 scan summation, KE=0.675eV, H2 ratio={:.3f}%'.format(cts[1])
plot_MS_bar(masses, cts, sttl, ttl)

#all data sets
dsAll = [ds352, ds349, ds356, ds349h2o]
dsp = [[ds[s]['sttl'], ds[s]['h2cts'].mean(), ds[s]['wcts'].mean() ]
		for ds in dsAll for s in sorted(ds.keys())]


#=================Consolidating H2/H2O leakag data============================#
# gather all datasets into one list
# iterate the list and

#H2 Thermal leakage
thermalH2Data = [ds349, ds352, ds356, ds018, ds028]

#=====================plotting SID Data for R1003============================#
dpath = inmsDataDir+'2015-272'
msSID = get_INMS_data(dpath, '*_MS_*')
mscSID = get_INMS_scan_data(msSID, id0='mass')

BKG = [0]
SID1 = range(1,8,2)
THM1 = range(2,9,2)
SID2 = range(9,16,2)
THM2 = range(10,17,2)
SID3 = range(17,25,2)
THM3 = range(18,25,2)

sid272 = {'SID1':SID1, 'THM1':THM1,
            'SID2':SID2, 'THM2':THM2,
            'SID3':SID3, 'THM3':THM3,
            'BKG':BKG}

ds272 = get_datasets(sid272)
sum_sets(ds272, mscSID)
dkeys = sorted(ds272.keys())

for k in sid272.keys():
    sid272[k] = sum(array([ds272[dkeys[i]]['sum'] for i in sid272[k]]),0)

M = mscSID[0][-1]
for k, s in sid272.iteritems():
    plot_MS_bar(M, s, k, 'R1003')

SID_1 = sid272['SID1']-sid272['THM1']
SID_2 = sid272['SID2']-sid272['THM2']
SID_3 = sid272['SID3']-sid272['THM3']

#=====================Plotting Simulted Energy Scan Data from GM=====================#
os.chdir(inmsDataDir)
simfile = glob.glob('*.dat')[0]
with open(simfile) as df:
        dlist = [item.strip('\r\n').split(',')
        			for item in df.readlines()]
simQL = array([d[0] for d in dlist], float)
simCnts = array([d[1] for d in dlist], float)
figure()
plot(simQL, simCnts/simCnts.max())
ylim(0, 1.1)
ylabel('Counts')
xlabel('QL3 Voltage')
title('Simulated Energy Scan 11.6eV CO2', fontsize=16)
draw()
savefig('Simulated Energy Scan 11.6eV CO2.png')
close()
        

#------------------------------------------------------------------------#
# Title: 	INMSAngularResponseDA.py
# Purpose:  Provide tools for reading, analyzing, and plotting INMS data.
#
#------------------------------------------------------------------------#
#
# How to use this software:
# Must have Python v2.7.5 or later and IPython installed and running.
# Open the IPython interpreter - QtConsole and run the following commands
#
# cd /Users/jgrimes/Documents/Ion Gun/data or relevant user directory
# from INMS_REU_Data_Analysis import *
# ion()
#
# All required scripts and packages have now been loaded and the directory
# has been set to the data location to access files. To run methods from
# the interpreter prompt type
# ...: functionName(arg1, arg2, ..., argN)
#------------------------------------------------------------------------#

#begin import packages
#from pylab import * #all of numpy and matplotlib
from numpy import *
from scipy import *
from datetime import datetime
from datetime import timedelta
from matplotlib.pyplot import *
from matplotlib.dates import DateFormatter
from scipy.interpolate import griddata

import os
import csv
import json
import glob

inmsDataDir = '/Users/jgrimes/Documents/Cassini/INMS Data/'
formatter = DateFormatter('%H:%M')

sfld = ['run', 'dataset', 'scan', 'mode',
        'initial_time', 'final_time']

mfld = ['ip','cs_fil1_emis_iset','cs_fil1_emis_imon', 'cs_anode1_imon',
        'cs_ionrep_imon', 'os_fil1_emis_imon', 'chamber_ion_pmon']

def get_INMS_data(datapath, srchstr):

    os.chdir(datapath)
    #make a list of appropriate files
    esFiles = glob.glob(srchstr)

    #read each file in the list
    esData = {}
    for dfile in esFiles:
        esData[dfile] = read_INMS_scan_file(dfile)

    return esData

def get_INMS_scan_data(data, id0='lens_vset', sfields='', mfields=''):
#
    if sfields == '':
        sfields = ['run', 'dataset', 'scan',
                    'initial_time', 'final_time']
    if mfields == '':
        mfields = ['cs_fil1_emis_imon', 'cs_anode1_imon',
                     'cs_ionrep_imon', 'chamber_ion_pmon']
    scanData = []
    for k in sorted(data.keys()):
        try:
            scan = data[k]['Scan']
            oneScan = [scan[item] for item in sfields]
            monitors = data[k]['Monitors']
            oneScan.extend([monitors[item] for item in mfields])
            titl = data[k]['Data'][0]
            cnts = data[k]['Data'][1][:, titl.index('counter1')]
            d = data[k]['Data'][1][:, titl.index(id0)]
            oneScan.extend([cnts, d])
            scanData.append(oneScan)
        except TypeError:
            print(k)
            pass
    return scanData

#===============DA tools for INMS Mass Scan Data============================#

def get_MS_scan_times(jn_file):
    tBkg, tSID, tTherm, SIDgrps = [],[],[],[]
    dset = 1
    with open(jn_file) as jn:
        for line in jn:
            if 'Scan' in line:
                jln0, jln1 = line.strip('[').split(']')
                if 'Background' in jln1:
                    tBkg.append(jln0)
                    SIDgrps.append({'Bkg':tBkg})
                elif 'SID' in jln1:
                    tSID.append(jln0)
                elif 'Thermal' in jln1:
                    tTherm.append(jln0)
            elif 'SID measurement completed' in line:
                tEnd = line.strip('[').split(']')[0]
                SIDgrps.append({'SID':tSID, 'Therm':tTherm, 'End':tEnd})
                tSID, tTherm = [], []
    return SIDgrps

def get_datasets(msFiles):
    dsets = {}
    for msf in msFiles:
        _,_,d,s = msf.strip('.txt').split('_')[2:]
        try:
            dsets[d].append(s)
        except KeyError:
            dsets[d] = [s]
    return dsets


def scale_ms_data(data, aveP=0, avEmis=0):
    if aveP==0:
        #if no average pressure provided, calculate for this set
        aveP = array([d[8] for d in data], dtype=float).mean()

    if avEmis==0:
        #if no average emission provided, calculate for this set
        avEmis = array([d[5] for d in data], dtype=float).mean()

    numscans = len(data)
    press = array([d[8] for d in data], dtype=float)
    normP = aveP/press
    emis = array([d[5] for d in data], dtype=float)
    normE = avEmis/emis
    scans = array([d[-2] for d in data])
    scaled = [scans[i]*normP[i]*normE[i] for i in range(0, numscans)]
    return scaled

def plot_MS_bar(M, D, sttl, ttl, yl='Counts', save=False, inline=True):
    figure()
    bar(M-.5, D)
    suptitle(sttl, fontsize=16)
    gca().set_title(ttl, fontsize=12)
    xlabel('Mass per Charge')
    ylabel(yl)
    xlim(.5, 25.5)
    ylim(0, 1.05*max(D))
    draw()
    savefig(sttl+'.png')

def plot_crosstalk(dsets, mscans):
    sum_sets(dsets, mscans)
    for s in dsets:
        plot_MS_bar(s['masses'], s['sum'], s['sttl'], s['ttl'])


def get_ds_ranges(dst):
    for k in dst.keys():
        sc = sorted(dst[k])
        i, f = int(sc[0][1:]), int(sc[-1][1:])
        dst[k] = {'srange':(i,f)}


def sum_sets(dsets, mscans):
    # sum the mass scans to create one composite mass scan for the dataset
    # extract all items relevant to scale and compare mass scan data
    for s in sorted(dsets.keys()):
        i, f = dsets[s]['srange']
        scanSum = sum(array([ms[-2] for ms in mscans[i-1: f]]),0)
        h2cts = array([ms[-2][1] for ms in mscans[i-1:f]])
        wcts = array([ms[-2][14] for ms in mscans[i-1:f]])
        csEmis = array([ms[8] for ms in mscans[i-1:f]])
        osEmis = array([ms[11] for ms in mscans[i-1:f]])
        chPress = array([ms[12] for ms in mscans[i-1:f]])
        run, mode, ip = mscans[f-1][0:7:3]
        masses = mscans[f-1][-1]
        numscans = f-(i-1)

        sttl = '{} - R{}'
        if 'bkg' in s:
            sttl =  '{} Background - R{}'
        elif mode=='OSNB':
            if float(mscans[f-1][7]) >= 1:
                sttl = '{} - CS On - R{}'
            else:
                sttl = '{} - CS Off - R{}'
        sttl = sttl.format(mode, run)
        if 'v_comp' in dsets[s].keys():
            sttl = sttl + ' {}kms'.format(dsets[s]['v_comp'])

        ttl = 'IP={}ms, {} scan summation, H2counts={:.0f}'
        ttl = ttl.format(ip, numscans, scanSum[1])

        dsets[s]['sum']=scanSum
        dsets[s]['h2cts']=h2cts
        dsets[s]['wcts'] = wcts
        dsets[s]['masses'] = masses
        dsets[s]['cemis'] = csEmis
        dsets[s]['oemis'] = osEmis
        dsets[s]['chPress'] = chPress
        dsets[s]['sttl']=sttl
        dsets[s]['ttl']=ttl


def compute_stats(dsets):
    pfmat = '{:<28}: sum={:12.0f}, std={:12.2f}, var={:12.2f}'
    hfmat = '{:<28}: {:>12}{:>12}{:>12}'
    print(hfmat.format('Dataset','H2 Sum','Avg. H2','Std. Dev.', 'Variance'))
    for dk in sorted(dsets.keys()):
        if 'h2cts' in dsets[dk].keys():
            ds = dsets[dk]
            ds['h2sum'] = ds['h2cts'].sum()
            ds['mean'] = ds['h2cts'].mean()
            ds['std'] = ds['h2cts'].std()
            ds['var'] = ds['h2cts'].var()
            pp = pfmat.format(ds['sttl'],ds['h2sum'],ds['mean'],ds['std'],ds['var'])
            print(pp)


def compute_h2ratios(dsets):
    #dsets must be a specific grouping of 5 data sets including, in order,
    #CSN background, OSNB background, OSNB-filOn, CSN, and OSNB-filOff.
    if len(dsets)==5:
        skeys = sorted(dsets.keys())

        #asign the variables
        bkgc, bkgo, osnb1, csn, osnb2 = [dsets[sk] for sk in skeys]

        #subtract the background
        csnb = csn['sum'] - bkgc['sum']
        osnb1b = osnb1['sum'] - bkgo['sum']
        osnb2b = osnb2['sum'] - bkgo['sum']

        #compute the ratios for all species
        osnb1r = osnb1b/csnb
        osnb2r = osnb2b/csnb

        #compute the uncertainty
        m1, m2, mc = osnb1['mean'], osnb2['mean'], csn['mean']
        m1r = m1/mc
        m2r = m2/mc
        sig1 = (m1/mc)*sqrt((1/m1)+(1/mc))
        sig2 = (m2/mc)*sqrt((1/m2)+(1/mc))

        #attach results to appropriate data sets
        osnb1['ratio'] = m1r
        osnb1['sigma'] = sig1

        osnb2['ratio'] = m2r
        osnb2['sigma'] = sig2

        #print results
        pfmat = '\tCS{}: \tR(H2)={:.2e}, uncertainty={:.2e}'
        print('\tH2 ratio OSNB/{}'.format(csn['sttl']))
        print(pfmat.format('On', m1r, sig1))
        print(pfmat.format('Off', m2r, sig2))

    else:
        print('Compute_h2ratios requires exactly 5 data sets.')


def compute_water_ratios(dsets):
    skeys = sorted(dsets.keys())
    mh1, mhc, mh2 = [dsets[sk]['mean'] for sk in skeys[2:]]
    mw1, mwc, mw2 = [dsets[sk]['sum'][14]/40 for sk in skeys[2:]]
    mhw1, mhw2 = mh1/mw1, mh2/mw2
    sig1 = (mh1/mw1)*sqrt((1/mh1)+(1/mw1))
    sig2 = (mh2/mw2)*sqrt((1/mh2)+(1/mw2))
    sig3 = mhc/mwc*sqrt((1/mhc)+(1/mwc))
    sig4 = mh1/mwc*sqrt((1/mh1)+(1/mhc))
    sig5 = mh2/mwc*sqrt((1/mh2)+(1/mhc))
    pfmat = '\t{}: \tR({})={:.2e}, uncertainty={:.2e}'
    print('\tH2/H2O ratio OSNB v_comp = 8.506 km/s')
    print(pfmat.format('OSNB CS-On', 'H2/H2O', mhw1, sig1))
    print(pfmat.format('OSNB CS-Off', 'H2/H2O',mhw2, sig2))
    print(pfmat.format('CS Only', 'H2/H2O', mhc/mwc, sig3))
    print(pfmat.format('OSNB CS_On', 'H2(OSNB)/H2O(CSN)', mh1/mwc, sig4))
    print(pfmat.format('OSNB CS_Off', 'H2(OSNB)/H2O(CSN)', mh2/mwc, sig5))

sfld = ['run', 'dataset', 'scan', 'mode',
        'initial_time', 'final_time']

mfld = ['ip','cs_fil1_emis_iset','cs_fil1_emis_imon', 'cs_anode1_imon',
        'cs_ionrep_imon', 'os_fil1_emis_imon', 'chamber_ion_pmon']


def get_crosstalk_data(dpath):
    msd = get_INMS_data(dpath, '*_MS_*')
    msc = get_INMS_scan_data(msd, id0='mass', sfields=sfld, mfields=mfld)
    dst = get_datasets(msd)
    get_ds_ranges(dst)
    sum_sets(dst, msc)
    compute_stats(dst)

    return dst, msd, msc


#sttl = 'CSN Background - R{}'
#ttlc = 'IP=1000ms, 40 scan summation, H2counts={:.0f}'.format(cts[1])

#sttl = 'OSNB Background - R{}'
#ttlc = 'IP=1000ms, 40 scan summation, H2counts={:.0f}'.format(cts[1])

#sttl = 'OSNB - CS On - R{}'
#ttlc = 'IP=1000ms, 40 scan summation, H2counts={:.0f}'.format(cts[1])

#sttl = 'CSN - R{}'
#ttlc = 'IP=1000ms, 40 scan summation, H2counts={:.0f}'.format(cts[1])

#sttl = 'OSNB - CS Off - R{}'
#ttlc = 'IP=1000ms, 40 scan summation, H2counts={:.0f}'.format(cts[1])

#sttl = 'Crosstalk Ratio, CS On - R{}'
#ttlr = 'IP=1000ms, 40 scan summation, H2 ratio = {:.3f}'.format(cts[1])

#sttl = 'Crosstalk Ratio, CS Off - R{}'
#ttlr = 'IP=1000ms, 40 scan summation, H2 ratio = {:.3f}'.format(cts[1])


#ct352 = [{dset:'bkgc', rrange}]