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