docs: def gp_ccX():
      """fit experimental data"""
      inDir, outDir = getWorkDirs()
      data, alldata = OrderedDict(), None
      for infile in os.listdir(inDir):
          # get key and import data
          key = os.path.splitext(infile)[0].replace('_', '/')
          data_import = np.loadtxt(open(os.path.join(inDir, infile), 'rb'))
          # convert to log10 vs log10 plot, z=log10(y) => dz=0.434*dy/y
          data_import[:,3] = 0.434*(data_import[:,3]/data_import[:,1])
          data_import[:,(0,1)] = np.log10(data_import[:,(0,1)])
          data_import[:,2] = 0
          # fill dictionary
          data[key] = data_import
          alldata = data[key] if alldata is None else np.vstack((alldata, data[key]))
      # fit linear part first
      lindata = alldata[alldata[:,0]>2.5]
      m = (lindata[-1,1]-lindata[0,1])/(lindata[-1,0]-lindata[0,0])
      t = lindata[0,1] - m * lindata[0,0]
      popt1, pcov = curve_fit(
          linear, lindata[:,0], lindata[:,1], p0=[m, t],
          sigma=lindata[:,3], absolute_sigma=True
      )
      # fit full range
      popt2, pcov = curve_fit(
          lambda x, c, d: fitfunc(x, popt1[0], popt1[1], c, d),
          alldata[:,0], alldata[:,1], sigma=alldata[:,3], absolute_sigma=True,
      )
      popt = np.hstack((popt1, popt2))
      model = lambda x: fitfunc(x, *popt)
      # calculate mean standard deviation of data from parameterization
      yfit = np.array([model(x) for x in alldata[:,0]])
      stddev = 1.5*np.sqrt( # multiple of "sigma"
          np.average((alldata[:,1]-yfit)**2, weights=1./alldata[:,3])
      )
      print 'stddev = %.2g' % stddev
      errorband = np.array([[x, model(x), 0, 0, stddev] for x in np.linspace(1,4)])
      # make plot
      fitdata = np.array([[x, model(x), 0, 0, 0] for x in np.linspace(1,4)])
      par_names = ['a', 'b', 'c', 'd']
      energies = [19.6, 27, 39, 62.4, 200]
      labels = dict(
          ('%s = %.3g' % (par_name, popt[i]), [3.3, 3-i*0.2, True])
          for i,par_name in enumerate(par_names)
      )
      ccX_vals = [10**model(np.log10(energy)) for energy in energies]
      ccX = [' '.join([
          '%g GeV:' % energy,
          '({})'.format(ufloat(ccX_vals[i], stddev/0.434*ccX_vals[i])),
          '{/Symbol \155}b'
      ]) for i,energy in enumerate(energies)]
      print ccX
      #labels.update(dict(
      #    (cc, [1+i*0.5, 4.5+(i%2+1)*0.2, True]) for i,cc in enumerate(ccX)
      #))
      make_plot(
          data = [errorband] + data.values() + [fitdata],
          properties = [
              'with filledcurves lt 1 lw 5 pt 0 lc %s' % default_colors[8]
          ] + [
              'lc %s lw 4 lt 1 pt 18 ps 1.5' % (default_colors[i])
              for i in xrange(len(data))
          ] + ['with lines lc 0 lw 4 lt 1'],
          titles = [''] + data.keys() + ['y = ax+b - e^{-cx+d}'],
          xlabel = 'x = log_{10}[{/Symbol \326}s_{NN} (GeV)]',
          ylabel = 'y = log_{10}[{/Symbol \163}@_{c@^{/=18-}c}^{NN} ({/Symbol \155}b)]',
          name = os.path.join(outDir, 'ccX'),
          size = '11.4in,8.3in', xr = [1, 4], yr = [0.5,4.5],
          key = ['bottom right', 'nobox', 'width -5'],
          tmargin = 0.98, rmargin = 0.99, bmargin = 0.13,
          lines = dict(
              ('y=%f' % (np.log10(energy)), 'lc 0 lw 2 lt 2') for energy in energies
          ), labels = labels,
      )