import warnings
import numpy as np
import time
from struct import unpack
from scipy.stats import norm
from scipy.signal import tukey
from Filter import median_filter
import Analysis, Filter, Constants

def savefits(data, filename, vmax=1.0, sps=1e6, bits=14, noise=False, clobber=True):
    """
    Save pulse/noise to FITS file
    """
    
    import pyfits as pf
    
    # Prepare data
    data = (np.asarray(data)/vmax*2**(bits-1)).round()
    
    # Column Name
    if noise:
    	colname = 'NoiseRec'
    else:
    	colname = 'PulseRec'
        
    # Columns
    col_t = pf.Column(name='TIME', format='1D', unit='s', array=np.zeros(data.shape[0], dtype=int))
    col_data = pf.Column(name=colname, format='%dI' % data.shape[1], unit='V', array=data)
    
    cols = pf.ColDefs([col_t, col_data])
    tbhdu = pf.BinTableHDU.from_columns(cols)
    
    # Name of extension
    exthdr = tbhdu.header
    exthdr['EXTNAME'] = ('Record', 'name of this binary table extension')
    exthdr['EXTVER'] = (1, 'extension version number')
    
    # Add more attributes
    exthdr['TSCAL2'] = (vmax/2**(bits-1), '[V/ch]')
    exthdr['TZERO2'] = (0., '[V]')
    exthdr['THSCL2'] = (sps**-1, '[s/bin] horizontal resolution of record')
    exthdr['THZER2'] = (0, '[s] horizontal offset of record')
    exthdr['THSAM2'] = (data.shape[1], 'sample number of record')
    exthdr['THUNI2'] = ('s', 'physical unit of sampling step of record')
    exthdr['TRMIN2'] = (-2**(bits-1)+1, '[channel] minimum number of each sample')
    exthdr['TRMAX2'] = (2**(bits-1)-1, '[channel] maximum number of each sample')
    exthdr['TRBIN2'] = (1, '[channel] default bin number of each sample')
    
    # More attributes
    exthdr['TSTART'] = (0, 'start time of experiment in total second')
    exthdr['TSTOP'] = (0, 'end time of experiment in total second')
    exthdr['TEND'] = (0, 'end time of experiment (obsolete)')
    exthdr['DATE'] = (time.strftime("%Y-%m-%dT%H:%M:%S", time.gmtime()), 'file creation date (UT)')
    
    # We anyway need Primary HDU
    hdu = pf.PrimaryHDU()
    
    # Write to FITS
    thdulist = pf.HDUList([hdu, tbhdu])
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        thdulist.writeto(filename, clobber=clobber)

def fopen(filename):
    """
    Read FITS file
    
    Parameters
    ==========
        filename:   file number to read
    
    Returns
    =======
        t:          time array
        wave:       waveform array
    """
    
    import pyfits as pf

    # Open fits file and get pulse/noise data
    header = pf.open(filename)
    wave = header[1].data.field(1).copy()
    dt = header[1].header['THSCL2']
    t = np.arange(wave.shape[-1]) * dt
    header.close()
    
    return t, wave

def yopen(filenumber, summary=False, nf=None, tmin=None, tmax=None, raw=False):
    """
    Read Yokogawa WVF file
    
    Parameters
    ==========
        filenumber: file number to read
        summary:    to summary waves (default: False)
        nf:         sigmas for valid data using median noise filter, None to disable noise filter (default: None)
        tmin:       lower boundary of time for partial extraction, scaler or list (Default: None)
        tmax:       upper boundary of time for partial extraction, scaler or list (Default: None)
        raw:        returns raw data without scaling/offsetting if True (Default: False)
    
    Returns
    =======
        if summary is False:
            [ t1, d1, t2, d2, t3, d3, ... ]
        
        if summary is True:
            [ t1, d1, err1, t2, d2, err2, ... ]
        
        if raw is True:
            t1 is a tuple of (hres1, hofs1, vres1, vofs1)
        
        where t1 is timing for 1st ch, d1 is data for 1st ch, err1 is error (1sigma) for 1st ch, and so on.
    """
    
    # Read header (HDR)
    h = open(str(filenumber) + ".HDR")
    lines = h.readlines()
    h.close()
    
    # Parse $PublicInfo
    for line in lines:
        token = line.split()
        
        if len(token) > 0:
            # Check endian
            if token[0] == "Endian":
                endian = '>' if token[1] == "Big" else '<'
            
            # Check data format
            if token[0] == "DataFormat":
                format = token[1]
                assert format == "Block"
                
            # Check # of groups
            if token[0] == "GroupNumber":
                groups = int(token[1])
            
            # Check # of total traces
            if token[0] == "TraceTotalNumber":
                ttraces = int(token[1])
            
            # Check data offset
            if token[0] == "DataOffset":
                offset = int(token[1])
    
    # Initialize containers
    traces = [None] * groups        # Number of traces for each group
    blocks = [None] * ttraces       # Number of blocks for each trace
    bsizes = [None] * ttraces       # Block size for each trace
    vres = [None] * ttraces         # VResolution for each trace
    voffset = [None] * ttraces      # VOffset for each trace
    hres = [None] * ttraces         # HResolution for each trace
    hoffset = [None] * ttraces      # HOffset for each trace
    
    # Parse $Group
    for line in lines:
        token = line.split()

        if len(token) > 0:
            # Read current group number
            if token[0][:6] == "$Group":
                cgn = int(token[0][6:]) - 1  # Current group number (minus 1)
            
            # Check # of traces in this group
            if token[0] == "TraceNumber":
                traces[cgn] = int(token[1])
                traceofs = np.sum(traces[:cgn], dtype=int)
                        
            # Check # of Blocks
            if token[0] == "BlockNumber":
                blocks[traceofs:traceofs+traces[cgn]] = [ int(token[1]) ] * traces[cgn]
            
            # Check Block Size
            if token[0] == "BlockSize":
                bsizes[traceofs:traceofs+traces[cgn]] = [ int(s) for s in token[1:] ]
            
            # Check VResolusion
            if token[0] == "VResolution":
                vres[traceofs:traceofs+traces[cgn]] = [ float(res) for res in token[1:] ]
            
            # Check VOffset
            if token[0] == "VOffset":
                voffset[traceofs:traceofs+traces[cgn]] = [ float(ofs) for ofs in token[1:] ]
            
            # Check VDataType
            if token[0] == "VDataType":
                assert token[1] == "IS2"
            
            # Check HResolution
            if token[0] == "HResolution":
                hres[traceofs:traceofs+traces[cgn]] = [ float(res) for res in token[1:] ]
            
            # Check HOffset
            if token[0] == "HOffset":
                hoffset[traceofs:traceofs+traces[cgn]] = [ float(ofs) for ofs in token[1:] ]
        
    # Data Initialization
    time = [ np.array(range(bsizes[t])) * hres[t] + hoffset[t] for t in range(ttraces) ]
    data = [ [None] * blocks[t] for t in range(ttraces) ]
    
    # Open WVF
    f = open(str(filenumber) + ".WVF", 'rb')
    f.seek(offset)
    
    # Read WVF
    if format == "Block":
        # Block format (assuming block size is the same for all the traces in Block format)
        for b in range(blocks[0]):
            for t in range(ttraces):
                if raw:
                    data[t][b] = np.array(unpack(endian + 'h'*bsizes[t], f.read(bsizes[t]*2)), dtype='int64')
                else:
                    data[t][b] = np.array(unpack(endian + 'h'*bsizes[t], f.read(bsizes[t]*2))) * vres[t] + voffset[t]
    else:
        # Trace format
        for t in range(ttraces):
            for b in range(blocks[t]):
                if raw:
                    data[t][b] = np.array(unpack(endian + 'h'*bsizes[t], f.read(bsizes[t]*2)), dtype='int64')
                else:
                    data[t][b] = np.array(unpack(endian + 'h'*bsizes[t], f.read(bsizes[t]*2))) * vres[t] + voffset[t]

    # Array conversion
    for t in range(ttraces):
        if raw:
            data[t] = np.array(data[t], dtype='int64')
        else:
            data[t] = np.array(data[t])
            
    
    # Tmin/Tmax filtering
    for t in range(ttraces):
        if type(tmin) == list or type(tmax) == list:
            if not (type(tmin) == list and type(tmax) == list and len(tmin) == len(tmax)):
                raise ValueError("tmin and tmax both have to be list and have to have the same length.")
            mask = np.add.reduce([ (time[t] >= _tmin) & (time[t] < _tmax) for (_tmax, _tmin) in zip(tmax, tmin)], dtype=bool)
        else:
            _tmin = np.min(time[t]) if tmin is None else tmin
            _tmax = np.max(time[t]) + 1 if tmax is None else tmax
            mask = (time[t] >= _tmin) & (time[t] < _tmax)
        
        data[t] = data[t][:, mask]
        time[t] = time[t][mask]
        
    f.close()
    
    if summary is False:
        # Return wave data as is
        if raw:
            return [ [ (hres[t], hoffset[t], vres[t], voffset[t]), data[t] ] for t in range(ttraces)  ]
        else:
            return [ [ time[t], data[t] ] for t in range(ttraces)  ]
    else:
        if nf is None:
            # Noise filter is off
            if raw:
                return [ [ (hres[t], hoffset[t], vres[t], voffset[t]), np.mean(data[t].astype(dtype='float64'), axis=0), np.std(data[t].astype(dtype='float64'), axis=0, ddof=1) ]
                            for t in range(ttraces) ]
            else:
                return [ [ time[t], np.mean(data[t], axis=0), np.std(data[t], axis=0, ddof=1) ]
                            for t in range(ttraces) ]
        else:
            # Noise filter is on
            if raw:
                return [ [ (hres[t], hoffset[t], vres[t], voffset[t]),
                            np.apply_along_axis(lambda a: np.mean(a[median_filter(a, nf)]), 0,  data[t].astype(dtype='float64')),
                            np.apply_along_axis(lambda a: np.std(a[median_filter(a, nf)], ddof=1), 0, data[t].astype(dtype='float64')) ]
                                for t in range(ttraces) ]
            else:
                return [ [ time[t],
                            np.apply_along_axis(lambda a: np.mean(a[median_filter(a, nf)]), 0,  data[t]),
                            np.apply_along_axis(lambda a: np.std(a[median_filter(a, nf)], ddof=1), 0, data[t]) ]
                                for t in range(ttraces) ]

def popen(filename, ch=None, raw=False):
    """
    Read pls file
    
    Parameters
    ==========
        filename:   file name to read
        ch:         returns data only for the given channel if given (Default: None)
        raw:        returns raw data without scaling/offsetting if True (Default: False)
    
    Returns
    =======
        if raw is True:
            [ header, vres, vofs, hres, hofs, tick, num, data, edata ]
        else:
            [ header, t, tick, num, data, edata ]
    """
    
    # Initialize
    header = {'COMMENT': []}
    vres = {}
    vofs = {}
    hres = {}
    hofs = {}
    tick = {}
    num  = {}
    data = {}
    edata = {}
    
    # Parser
    def parser():
        """
        PLS Data Parser (generator)
        """
        
        # Initialization
        samples = -1
        extra = 0
        chunk = ''
        isHeader = True
        
        while True:
            while len(chunk) < 2:
                chunk += yield
        
            # Get the magic character
            magic = chunk[0]

            if isHeader and magic == 'C':
                # Comment
                while len(chunk) < 80:
                    chunk += yield
                header['COMMENT'].append(chunk[2:80])
                chunk = chunk[80:]
            
            elif isHeader and magic == 'V':
                # Version
                while len(chunk) < 80:
                    chunk += yield
                header['VERSION'] = chunk[2:80]
                chunk = chunk[80:]
            
            elif isHeader and magic == 'O':
                # Date
                while len(chunk) < 10:
                    chunk += yield
                _m, _d, _y = map(int, chunk[2:10].split())
                header['DATE'] = "%d/%d/%d" % (_y, _m, _d)
                chunk = chunk[10:]

            elif isHeader and magic == 'S':
                # Number of Samples
                while len(chunk) < 7:
                    chunk += yield
                header['SAMPLES'] = samples = int(chunk[2:7])
                chunk = chunk[7:]

            elif isHeader and magic == 'E':
                # Extra Bytes
                while len(chunk) < 7:
                    chunk += yield
                header['EXTRA'] = extra = int(chunk[2:7])
                chunk = chunk[7:]
                
            elif isHeader and magic == 'P':
                # Discriminator
                while len(chunk) < 78:
                    chunk += yield
                _dis = chunk[2:78].split()
                if _dis[0] == '01':
                    header['ULD'] = eval(_dis[1])
                elif _dis[0] == '02':
                    header['LLD'] = eval(_dis[1])
                chunk = chunk[78:]
            
            elif isHeader and magic == 'N':
                # Normalization
                while len(chunk) < 47:
                    chunk += yield
                _ch, _hofs, _hres, _vofs, _vres = chunk[2:47].split()
                _ch = int(_ch)
                vres[_ch] = eval(_vres)
                vofs[_ch] = eval(_vofs)
                hres[_ch] = eval(_hres)
                hofs[_ch] = eval(_hofs)
                chunk = chunk[47:]
            
            elif magic == 'D':
                # Data
                isHeader = False
                
                if samples < 0:
                    raise ValueError("Invalid number of samples.")
                while len(chunk) < (11 + samples*2):
                    chunk += yield
                _ch, _tick, _num = unpack('<BII', chunk[2:11])
                
                if not data.has_key(_ch):
                    data[_ch] = bytearray()
                    tick[_ch] = []
                    num[_ch] = []
                    edata[_ch] = bytearray()
                
                data[_ch] += chunk[11:11 + samples*2]
                tick[_ch].append(_tick)
                num[_ch].append(_num)
                edata[_ch] += chunk[11 + samples*2:11 + samples*2 + extra]
                
                chunk = chunk[11 + samples*2 + extra:]
                
            else:
                # Skip unknown magic
                chunk = chunk[1:]
    
    # Open pls file and read by chunks
    f = open(filename, 'rb')

    # Start parser
    p = parser()
    p.next()
    
    # Read by chunk and parse it
    with open(filename, 'rb') as f:
        while True:
            chunk = f.read(1024*1024)   # read 1 MB
            if not chunk:
                break
            p.send(chunk)

    # Convert buffer to numpy array
    for k in ([ch] if ch else data.keys()):
        data[k] = np.frombuffer(data[k], dtype='>i2').reshape(-1, header['SAMPLES'])
        edata[k] = np.frombuffer(edata[k], dtype='>u1').reshape(-1, header['SAMPLES'])
    
    if raw:
        if ch:
            return header, vres[ch], vofs[ch], hres[ch], hofs[ch], tick[ch], num[ch], data[ch], edata[ch]
        else:
            return header, vres, vofs, hres, hofs, tick, num, data, edata
    else:
        t = {}

        for k in ([ch] if ch else data.keys()):
            # Normalize data using res/ofs
            t[k] = (np.arange(header['SAMPLES']) + hofs[k]) * hres[k]
            data[k] = (np.asarray(data[k]) + vofs[k]) * vres[k]
        
        if ch:
            return header, t[ch], tick[ch], num[ch], data[ch], edata[ch]
        else:
            return header, t, tick, num, data, edata

def tesana(t, p, n, lpfc=None, hpfc=None, binsize=1, max_shift=10,
            thre=0.4, filt=None, nulldc=False, offset=False, center=False, sigma=3,
            gain=None, dsr=None, shift=False, ocmethod="ols", flip=False, atom="Mn",
            kbfit=False, ignorekb=False, method="mle",
            rshunt=None, tbias=None, ites=None, ka_min=80, kb_min=40,
            tex=False, plotting=True, savedat=False, session="Unnamed"):
    """
    Perform TES Analysis
    
    Parameters (and their default values):
        t:          time data (array-like)
        p:          pulse data (array-like)
        n:          noise data (array-like)
        lpfc:       low-pass filter cut-off frequency in bins (Default: None)
        hpfc:       high-pass filter cut-off frequency in bins (Default: None)
        binsize:    energy bin size for histograms and fittings (only for ls ans cs) in eV (Default: 1)
        max_shift:  maximum allowed shifts to calculate maximum cross correlation (Default: 10)
        thre:       correlation threshold for offset correction (Default: 0.4)
        filt:       window function (hanning/hamming/blackman/tukey) (Default: None)
        nulldc:     nullify the DC bin when template generation (Default: False)
        offset:     subtract DC offset (Default: False)
        center:     centering pulse rise (Default: False)
        sigma:      sigmas for median filter (Default: 3)
        gain:       feedback gain for current-space conversion (Default: None)
        dsr:        down-sampling rate (Default: None)
        shift:      treat dE as energy shift instead of scaling (Default: False)
        ocmethod:   offset correction fitting method (ols/odr) (Default: ols)
        flip:       flip x and y when offset correction fitting (Default: False)
        atom:       atom to fit (Default: Mn)
        kbfit:      fit Kb line (Default: False)
        ignorekb:   ignore Kb line when linearity correction (Default: False)
        method:     fitting method (mle/ls/cs) (Default: mle)
        rshunt:     shunt resistance value for r-space conversion (Default: None)
        tbias:      TES bias current for r-space conversion (Default: None)
        ites:       TES current for r-space conversion (Default: None)
        ka_min:     minimum counts to group bins for Ka line (valid only for ls/cs fittings) (Default: 80)
        ka_min:     minimum counts to group bins for Kb line (valid only for ls/cs fittings) (Default: 20)
        tex:        use TeX for plots (Default: False)
        plotting:   generate and save plots (Default: True)
        savedat:    save data to files (Default: False)
        session:    session name for plots and data files (Default: Unnamed)
        
    Note:
        - Use offset option when using filt option
        - Consider using center option when using filt option
    """

    if plotting:
        # Import matplotlib
        import matplotlib
        matplotlib.use('Agg')
        matplotlib.rcParams['text.usetex'] = str(tex)
        from pylab import figure, plot, errorbar, hist, axvline, xlim, ylim, loglog, xlabel, ylabel, legend, tight_layout, savefig
        
        print "Session: %s" % session
    
    # Preparation
    p = np.asarray(p)
    n = np.asarray(n)
    t = np.asarray(t)
    dt = np.diff(t)[0]
    df = (dt * t.shape[-1])**-1
    
    # Subtract offset
    if offset:
        ofs = np.median(n)
        p -= ofs
        n -= ofs
    
    # Convert to current-space if needed
    if gain:
        print "Converting to current-space"
        p /= gain
        n /= gain
    
    # Convert to resistance-space
    Rspace = False
    if gain and rshunt and tbias and ites:
        print "Converting to resistance-space"
        ofs = np.median(n)
        p += (ites - ofs)
        n += (ites - ofs)
        
        # Convert to resistance
        p = (tbias - p) * rshunt / p
        n = (tbias - n) * rshunt / n
        
        Rspace = True
    
    # Down-sample
    if dsr > 1:
        p = p[:,:p.shape[-1]/dsr*dsr].reshape(p.shape[0], -1, dsr).mean(axis=-1)
        n = n[:,:n.shape[-1]/dsr*dsr].reshape(n.shape[0], -1, dsr).mean(axis=-1)
        dt *= dsr
        t = t[::dsr]

    # Pulse centering (for filtering)
    if center:
        # Roll pulse to the center
        r = p.shape[-1] / 2 - np.median(abs(p - Filter.offset(p)[:, np.newaxis]).argmax(axis=-1))
        p = np.hstack((p[...,-r:], p[...,:-r]))

    # Calculate offset (needs to be done before applying filter)
    if p.size > 0:
        offset = Filter.offset(p)
    
    # Generate Filter
    if filt is None:
        pass
    else:
        if filt.lower() == "hanning":
            f = np.hanning(p.shape[-1])
        elif filt.lower() == "hamming":
            f = np.hamming(p.shape[-1])
        elif filt.lower() == "blackman":
            f = np.blackman(p.shape[-1])
        elif filt.lower() == "tukey":
            f = tukey(p.shape[-1])
        else:
            raise ValueError('Unsupported filter: %s' % filt.lower())
            
        print "Window filter function: %s" % filt.lower()
        
        # Amplitude correction
        cf = f.sum() / len(f)
        p *= (f / cf)
        n *= (f / cf)

        # Equivalent noise bandwidth correction
        enb = len(f)*(f**2).sum()/f.sum()**2
        df *= enb
        
    if p.size > 0:
        # Calculate averaged pulse
        avgp = Filter.average_pulse(p, max_shift=max_shift)

        if savedat:
            np.savetxt('%s-averagepulse.dat' % session, np.vstack((t, avgp)).T,
                header='Time (s), Averaged Pulse (%s)' % ('R' if Rspace else ('A' if gain else 'V')), delimiter='\t')

        if plotting:
            figure()
            plot(t, avgp)
            xlabel('Time$\quad$(s)')
            ylabel('Averaged Pulse$\quad$(%s)' % ('R' if Rspace else ('A' if gain else 'V')))
            tight_layout()
            savefig('%s-averagepulse.pdf' % session)

        # Calculate averaged pulse spectrum
        avgps = np.sqrt(Filter.power(avgp)) / df

        if savedat:
            np.savetxt('%s-avgpulse-power.dat' % session, np.vstack((np.arange(len(avgps))*df, avgps)).T,
                header='Frequency (Hz), Average Pulse Power (%s/srHz)' % ('R' if Rspace else ('A' if gain else 'V')), delimiter='\t')
    
        if plotting:
            avgps[0] = 0    # for better plot
            figure()
            plot(np.arange(len(avgps))*df, avgps)
            loglog()
            xlabel('Frequency$\quad$(Hz)')
            ylabel('Average Pulse Power$\quad$(%s/Hz)' % ('R' if Rspace else ('A' if gain else 'V')))
            tight_layout()
            savefig('%s-avgpulse-power.pdf' % session)

    if n.size > 0:
        # Plot noise spectrum
        avgns = np.sqrt(Filter.average_noise(n) / df)

        if savedat:
            np.savetxt('%s-noise.dat' % session, np.vstack((np.arange(len(avgns))*df, avgns)).T,
                header='Frequency (Hz), Noise (%s/srHz)' % ('R' if Rspace else ('A' if gain else 'V')), delimiter='\t')
    
        if plotting:
            avgns[0] = 0    # for better plot
            figure()
            plot(np.arange(len(avgns))*df, avgns)
            loglog()
            xlabel('Frequency$\quad$(Hz)')
            ylabel('Noise$\quad$(%s/$\sqrt{\mathrm{Hz}}$)' % ('R' if Rspace else ('A' if gain else 'V')))
            tight_layout()
            savefig('%s-noise.pdf' % session)

    if p.size > 0 and n.size > 0:
        # Generate template
        tmpl, sn = Filter.generate_template(p, n, lpfc=lpfc, hpfc=hpfc, nulldc=nulldc, max_shift=max_shift)

        if savedat:
            np.savetxt('%s-template.dat' % session, np.vstack((t, tmpl)).T,
                header='Time (s), Template (A.U.)', delimiter='\t')
            np.savetxt('%s-sn.dat' % session, np.vstack((np.arange(len(sn))*df, sn/np.sqrt(df))).T,
                header='Frequency (Hz), S/N (/srHz)', delimiter='\t')
    
        if plotting:
            # Plot template
            figure()
            plot(t, tmpl)
            xlabel('Time$\quad$(s)')
            ylabel('Template$\quad$(A.U.)')
            tight_layout()
            savefig('%s-template.pdf' % session)
        
            # Plot SNR
            figure()
            plot(np.arange(len(sn))*df, sn/np.sqrt(df))
            loglog()
            xlabel('Frequency$\quad$(Hz)')
            ylabel('S/N$\quad$(/$\sqrt{\mathrm{Hz}}$)')
            tight_layout()
            savefig('%s-sn.pdf' % session)
    
        # Calculate baseline resolution
        print "Resolving power: %.2f (%.2f eV @ 5.9 keV)" % (np.sqrt((sn**2).sum()*2), Analysis.baseline(sn))
    
        # Perform optimal filtering
        pha_p = Filter.optimal_filter(p, tmpl, max_shift=max_shift)
        pha_n = Filter.optimal_filter(n, tmpl, max_shift=0)
    
        # Offset correction
        (a, b), coef = Analysis.fit_offset(pha_p, offset, sigma=sigma, method=ocmethod, flip=flip)
    
        if coef > thre:
            oc_pha_p = Analysis.offset_correction(pha_p, offset, b)
            oc_pha_n = Analysis.offset_correction(pha_n, offset, b)
            print "Offset correction with: PHA = %f * (1 + %f * Offset)" % (a, b)
    
            if plotting:
                figure()
                ka = Analysis.ka(np.vstack((pha_p, offset)).T, sigma=sigma)
                plot(ka.T[1], ka.T[0], '.', c='k')
                x_min, x_max = xlim()
                ofs = np.linspace(x_min, x_max)
                label = '$\mathrm{PHA}=%.2f\\times(1+%.2f\\times\mathrm{Offset})$' % (a, b)
                plot(ofs, a*(1+b*ofs), 'r-', label=label)
                xlabel('Offset$\quad$(V)')
                ylabel('PHA$\quad$(V)')
                legend(frameon=False)
                tight_layout()
                savefig('%s-offset.pdf' % session)
        else:
            oc_pha_p = pha_p
            oc_pha_n = pha_n
            print "Skipped offset correction: correlation coefficient (%f) is too small" % coef

        # Check line database
        if "%sKa" % atom not in Constants.LE.keys() or "%sKb" % atom not in Constants.LE.keys():
            raise ValueError('Unsupported atom: %s' % atom)

        # Linearity correction
        pha_line_center = np.asarray([ np.median(Analysis.ka(oc_pha_p, sigma=sigma)), np.median(Analysis.kb(oc_pha_p, sigma=sigma)) ])
        line_energy = np.asarray([ Constants.LE['%sKa' % atom], Constants.LE['%sKb' % atom] ])

        if ignorekb:
            a, b = Analysis.fit_linearity([pha_line_center[0]], [line_energy[0]], deg=1)
            print "Linearity correction with: PHA = %e * E" % (b)
        else:
            a, b = Analysis.fit_linearity(pha_line_center, line_energy, deg=2)
            print "Linearity correction with: PHA = %e * E^2 + %e * E" % (a, b)
            print "MnKb saturation ratio: %.2f %%" % ((pha_line_center[1]/pha_line_center[0])/(line_energy[1]/line_energy[0])*100)

        lc_pha_p = Analysis.linearity_correction(oc_pha_p, a, b)
        lc_pha_n = Analysis.linearity_correction(oc_pha_n, a, b)

        if savedat:
            np.savetxt('%s-linearity.dat' % session, array([pha_line_center[0]]) if ignorekb else pha_line_center[np.newaxis,:],
                header='%sKa PHA' % atom if ignorekb else '%sKa PHA, %sKb PHA' % (atom, atom), delimiter='\t')
    
        if plotting:
            figure()
            x = np.linspace(0, 7e3)
            if ignorekb:
                plot(line_energy[0]/1e3, pha_line_center[0], '+', color='b')
                plot(x/1e3, x*b, 'r--')
            else:
                plot(line_energy/1e3, pha_line_center, '+', color='b')
                plot(x/1e3, x**2*a+x*b, 'r--')
            xlim((0, 7))
            xlabel('Energy$\quad$(keV)')
            ylabel('PHA$\quad$(a.u.)')
            tight_layout()
            savefig('%s-linearity.pdf' % session)
    
        # Energy Spectrum
        if plotting:
            figure()
            hcount, hbin, hpatch = hist(lc_pha_p[lc_pha_p==lc_pha_p]/1e3, bins=7000/binsize, histtype='stepfilled', color='y')
            xlim(0, 7)
            xlabel('Energy$\quad$(keV)')
            ylabel('Count')
            tight_layout()
            savefig('%s-spec.pdf' % session)

        if savedat:
            hcount, hbin = np.histogram(lc_pha_p[lc_pha_p==lc_pha_p]/1e3, bins=7000/binsize)
            np.savetxt('%s-spec.dat' % session, np.vstack(((hbin[1:]+hbin[:-1])/2, hcount)).T,
                header='Energy (keV), Count', delimiter='\t')
    
        # Line fitting
        def _line_fit(data, min, line):

            # Fit
            (dE, width), (dE_error, width_error), e = Analysis.fit(data, binsize=binsize, min=min, line=line, shift=shift, method=method)

            if method == "cs":
                chi_squared, dof = e

            if method in ("mle", "ls"):
                print "%s: %.2f +/- %.2f eV @ Ec%+.2f eV" \
                    % (line, width, width_error, dE)
            elif method == "cs":
                print "%s: %.2f +/- %.2f eV @ Ec%+.2f eV (Red. chi^2 = %.1f/%d = %.2f)" \
                    % (line, width, width_error, dE, chi_squared, dof, chi_squared/dof)

            return dE, width, width_error
    
        def _line_spectrum(data, min, line, dE, width, width_error):

            # Draw histogram
            n, bins = Analysis.histogram(data, binsize=binsize)

            if method in ("cs"):
                gn, gbins = Analysis.group_bin(n, bins, min=min)
            else:
                # No grouping in mle and ls
                gn, gbins = n, bins

            ngn = gn/(np.diff(gbins))
            ngn_sigma = np.sqrt(gn)/(np.diff(gbins))
            cbins = (gbins[1:]+gbins[:-1])/2

            if plotting:
                figure()

                if width_error is not None:
                    label = 'FWHM$=%.2f\pm %.2f$ eV' % (width, width_error)
                else:
                    label = 'FWHM$=%.2f$ eV (Fixed)' % width

                if method == "cs":
                    errorbar(cbins, ngn, yerr=ngn_sigma, xerr=np.diff(gbins)/2, capsize=0, ecolor='k', fmt=None, label=label)
                else:
                    hist(data, bins=gbins, weights=np.ones(len(data))/binsize, histtype='step', ec='k', label=label)

                E = np.linspace(bins.min(), bins.max(), 1000)

                model = Analysis.normalization(ngn, gbins, dE, width, line=line, shift=shift) \
                            * Analysis.line_model(E, dE, width, line=line, shift=shift, full=True)

                # Plot theoretical model
                plot(E, model[0], 'r-')

                # Plot fine structures
                for m in model[1:]:
                    plot(E, m, 'b--')

                xlabel('Energy$\quad$(eV)')
                ylabel('Normalized Count$\quad$(count/eV)')
                legend(frameon=False)

                ymin, ymax = ylim()
                ylim(ymin, ymax*1.1)
                tight_layout()

                savefig("%s-%s.pdf" % (session, line))

            if savedat:
                np.savetxt('%s-%s.dat' % (session, line), np.vstack((cbins, gn)).T,
                    header='Energy (keV), Count', delimiter='\t')
    
        ## Ka
        ka = Analysis.ka(lc_pha_p, sigma=sigma)
        dE, width, width_error = _line_fit(ka, ka_min, "%sKa" % atom)
        _line_spectrum(ka, ka_min, "%sKa" % atom, dE, width, width_error)

        ## Kb
        kb = Analysis.kb(lc_pha_p, sigma=sigma)
        if kbfit:
            dE, width, width_error = _line_fit(kb, kb_min, "%sKb" % atom)
        else:
            width_error = None
        _line_spectrum(kb, kb_min, "%sKb" % atom, dE, width, width_error)

        ## Baseline
        f_pha_n = lc_pha_n[Filter.median_filter(lc_pha_n, sigma=sigma)]
        baseline = Analysis.sigma2fwhm(np.std(f_pha_n))
        print "Baseline resolution: %.2f eV" % baseline

        n, bins = Analysis.histogram(f_pha_n, binsize=binsize)
    
        if savedat:
            np.savetxt('%s-baseline.dat' % session, np.vstack(((bins[1:]+bins[:-1])/2, n)).T,
                header='Energy (keV), Count', delimiter='\t')

        if plotting:
            figure()
            label = 'FWHM$=%.2f$ eV' % baseline
            hist(f_pha_n, bins=bins, weights=np.ones(len(f_pha_n))/binsize, histtype='step', ec='k', label=label)
            mu, sigma = norm.fit(f_pha_n)
            E = np.linspace(bins.min(), bins.max(), 1000)
            plot(E, norm.pdf(E, loc=mu, scale=sigma)*len(f_pha_n), 'r-')
        
            xlabel('Energy$\quad$(eV)')
            ylabel('Normalized Count$\quad$(count/eV)')
        
            legend(frameon=False)
        
            tight_layout()
        
            savefig('%s-baseline.pdf' % session)