docs: def gp_rdiff(version, nomed, noxerr, diffRel, divdNdy):
    """example for ratio or difference plots using QM12 data (see gp_panel)
  
    - uses uncertainties package for easier error propagation and rebinning
    - stat. error for medium = 0!
    - stat. error for cocktail ~ 0!
    - statistical error bar on data stays the same for diff
    - TODO: implement ratio!
    - TODO: adjust statistical error on data for ratio!
    - TODO: adjust name and ylabel for ratio
  
    .. image:: pics/diffAbsQM12.png
       :width: 450 px
  
    :param version: plot version
    :type version: str
    :param nomed: don't plot medium
    :type nomed: bool
    :param noxerr: don't plot x-errors
    :type noxerr: bool
    """
    inDir, outDir = getWorkDirs()
    inDir = os.path.join(inDir, version)
    data, cocktail, medium, rhofo, vacrho = \
            OrderedDict(), OrderedDict(), OrderedDict(), OrderedDict(), OrderedDict()
    #scale = { # QM14 (19 GeV skip later, factor here only informational)
    #  '19.6': 1.0340571932983775, '200': 1.0, '39': 0.7776679085382481,
    #  '27': 0.6412140408244136, '62.4': 0.9174700031778402
    #}
    scale = {
        '19.6': 1.3410566491548412, '200': 1.1051002240771077,
        '39': 1.2719203877292842, '27': 1.350873678084769,
        '62.4': 1.2664666321635087
    }
    yunit = 1.0e-3 if not diffRel else 1.
    for infile in os.listdir(inDir):
      if infile == "cocktail_contribs": continue
      energy = re.compile('\d+').search(infile).group()
      data_type = re.sub('%s\.dat' % energy, '', infile)
      energy = getEnergy4Key(energy)
      file_url = os.path.join(inDir, infile)
      data_import = np.loadtxt(open(file_url, 'rb'))
      if data_type != 'data' and (
          (version == 'QM14' and energy != '19.6') or version == 'LatestPatrickJieYi'
      ):
         data_import[:,(1,3,4)] /= scale[energy]
      if version == 'LatestPatrickJieYi':
          if data_type == 'data':
              data_import = data_import[(data_import[:,0] > 0.14) & (data_import[:,0] < 1.0)]
          else:
              data_import = data_import[data_import[:,0] < 1.0]
      if data_type == 'data': data[energy] = data_import
      elif data_type == 'cocktail': cocktail[energy] = data_import
      elif data_type == 'rho' or data_type == 'vacRho' or data_type == 'medium':
          if noxerr and not diffRel: data_import[:,2:] = 0.
          data_import[:,1] /= yunit
          if data_type == 'rho':
              mask = data_import[:,1] > 0.1
              rhofo[energy] = data_import if diffRel else data_import[mask]
          elif data_type == 'vacRho':
              mask = (data_import[:,0] > 0.35) & (data_import[:,1] > 0.01)
              vacrho[energy] = data_import if diffRel else data_import[mask]
          elif not nomed and data_type == 'medium':
              medium[energy] = data_import
    nSetsData = len(data)
  
    shift = { '19.6': '1e0', '27': '1e1', '39': '1e2', '62.4': '1e3', '200': '1e4'
    } if not diffRel else {
        '19.6': '1', '27': '8', '39': '50', '62.4': '200', '200': '900'
    }
    dataOrdered = OrderedDict()
    for energy in sorted(data, key=float, reverse=True):
      # data & bin edges
      # getUArray propagates stat/syst errors separately internally but
      # errors need to be doubled to retrieve correct errors
      uData = getUArray(data[energy])
      eData = getEdges(data[energy])
      uCocktail = getUArray(cocktail[energy])
      eCocktail = getEdges(cocktail[energy])
      loop = [eData]
      if energy in medium and diffRel:
        uMedium = getUArray(medium[energy])
        eMedium = getEdges(medium[energy])
        loop.append(eMedium)
      if energy in rhofo and diffRel:
        uRho = getUArray(rhofo[energy])
        eRho = getEdges(rhofo[energy])
        loop.append(eRho)
      if energy in vacrho and diffRel:
        uVacRho = getUArray(vacrho[energy])
        eVacRho = getEdges(vacrho[energy])
        loop.append(eVacRho)
      # loop data/medium bins
      for l, eArr in enumerate(loop):
        for i, (e0, e1) in enumzipEdges(eArr):
          logging.debug('%s/%d> %g - %g:' % (energy, l, e0, e1))
          # get cocktail sum in data bin range
          # value+/-0.5*tot.uncert.
          uCocktailSum = getCocktailSum(e0, e1, eCocktail, uCocktail)
          if uCocktailSum == 0.: continue
          # calc. difference and divide by data binwidth again
          # + set data point
          if l == 0:
            uDiff = uData[i] # value+/-0.5*tot.uncert.
            if diffRel:
              uDiff /= uCocktailSum # value+/-0.5*tot.uncert.
            else:
              uDiff -= uCocktailSum
              uDiff /= data[energy][i,2] * 2 * yunit
            dp = [
              data[energy][i,0], uDiff.nominal_value,
              data[energy][i,2] if not noxerr else 0.,
              getErrorComponent(uDiff, 'stat'),
              getErrorComponent(uDiff, 'syst')
            ]
            key = ' '.join([energy, 'GeV'])
            if noxerr:
                if diffRel:
                    key += ' {/Symbol \264} %s' % shift[energy]
                else:
                    expon = shift[energy].split('e')[1]
                    key += ' {/Symbol \264} 10^{%s}' % expon
          elif l == 1:
            # only done if diffRel
            uDiff = uMedium[i]
            uDiff /= uCocktailSum
            dp = [
              medium[energy][i,0], uDiff.nominal_value+1,
              medium[energy][i,2] if not noxerr else 0.,
              0., 0. # both errors included in data points
            ]
            key = ' '.join([energy, 'GeV (Med.)'])
          elif l == 2:
            # only done if diffRel
            uDiff = uRho[i]
            uDiff /= uCocktailSum
            dp = [
              rhofo[energy][i,0], uDiff.nominal_value+1.,
              rhofo[energy][i,2] if not noxerr else 0.,
              0., 0. # both errors included in data points
            ]
            key = ' '.join([energy, 'GeV (RhoFO.)'])
          elif l == 3:
            # only done if diffRel
            uDiff = uVacRho[i]
            uDiff /= uCocktailSum
            dp = [
              vacrho[energy][i,0], uDiff.nominal_value+1.,
              vacrho[energy][i,2] if not noxerr else 0.,
              0., 0. # both errors included in data points
            ]
            key = ' '.join([energy, 'GeV (VacRho.)'])
          # build list of data points
          if diffRel or l == 0:
            if dp[0] > 0.7425 and dp[0] < 0.825: continue # mask out omega region
            if dp[0] > 0.97 and dp[0] < 1.0495: continue # mask out phi region
          if key in dataOrdered:
              dataOrdered[key] = np.vstack([dataOrdered[key], dp])
          else:
              dataOrdered[key] = np.array([ dp ])
      if not diffRel:
        if energy in medium:
          dataOrdered[' '.join([energy, 'GeV (Med.)'])] = medium[energy]
        if energy in rhofo:
          dataOrdered[' '.join([energy, 'GeV (RhoFO.)'])] = rhofo[energy]
        if energy in vacrho:
          dataOrdered[' '.join([energy, 'GeV (VacRho.)'])] = vacrho[energy]
  
    # make plot
    nSets = len(dataOrdered)
    nCats = 4
    nSetsPlot = nSets/nCats if nSets > nSetsData else nSets
    props = [
      'lt 1 lw 4 ps 1.5 lc %s pt 18' % default_colors[i]
        for i in reversed(range(nSetsPlot))
    ]
    titles = dataOrdered.keys()
    if nSets > nSetsData:
      props = zip_flat(props, *[
          [
              'with lines lt %d lw 4 lc %s' % (j+1, default_colors[i])
              for i in reversed(range(nSetsPlot))
          ]
          for j in xrange(nCats-1)
      ])
      titles = zip_flat(dataOrdered.keys()[::nCats], *[ [''] * nSetsPlot for j in xrange(nCats-1) ])
    global labels
    labels = {
      '{BES: STAR Preliminary}' if version == 'QM12Latest200' or \
        version == 'QM14' or version == 'LatestPatrickJieYi'
      else 'STAR Preliminary': [
          0.4 if diffRel else 0.2,0.09 if not diffRel and noxerr else 0.75,False
      ],
      '{200 GeV: PRL 113 022301' if version == 'QM12Latest200' \
        or version == 'QM14' or version == 'LatestPatrickJieYi'
      else '': [0.4 if diffRel else 0.2,0.04 if not diffRel and noxerr else 0.7,False],
    }
    yr = [.6,2.5e3] if diffRel else [0.05,1.5e5]
    if noxerr:
        for k,d in dataOrdered.iteritems():
            energy = getEnergy4Key(re.compile('\d+').search(k).group())
            d[:,(1,3,4)] *= float(shift[energy])
    gpcalls = [
        'object 1 rectangle back fc rgb "grey" from 0.75,%f to 0.825,%f' % \
        (1.7 if diffRel else 0.5, yr[1]),
        'object 2 rectangle back fc rgb "grey" from 0.96,%f to 1.0495,%f' % \
        (1.7 if diffRel else 0.5, yr[1]),
        'object 3 rectangle back fc rgb "#C6E2FF" from 0.4,%f to 0.75,%f' % \
        (1.7 if diffRel else 0.5, yr[1]),
        'boxwidth 0.01 absolute',
    ]
    hline = 1. if diffRel else .5
    lines = dict(
        (('x=%g' % (hline*float(shift[energy]))), 'lc rgb "black" lw 4 lt 4')
        for energy in shift
    )
    pseudo_point = np.array([[-1,1,0,0,0]])
    make_plot(
      data = dataOrdered.values() + [
          pseudo_point, pseudo_point, pseudo_point, pseudo_point
      ],
      properties = props + [
          'with lines lt %d lw 4 lc rgb "black"' % (lt+1)
          for lt in xrange(nCats)
      ],
      titles = titles + [
          'HMBT + QGP', 'BW/FO-{/Symbol \162}', '{/Symbol \162}/{/Symbol \167} VacSF+FB+FO',
          'baseline', #'%g%s' % (hline, ' {/Symbol \264} 10^{-3}' if not diffRel else '')
      ],
      name = os.path.join(outDir, 'diff%s%s%s%s' % (
        'Rel' if diffRel else 'Abs', version,
        'NoMed' if nomed else '', 'NoXErr' if noxerr else ''
      )),
      xlabel = 'dielectron invariant mass, M_{ee} (GeV/c^{2})',
      ylabel = 'Enhancement Ratio' if diffRel else 'Excess Yield / dM_{ee} ({/Symbol \264} 10^{-3} (GeV/c^2)^{-1})',
      #labels = labels,
      xr = [0.18,0.97], yr = yr, ylog = True,
      key = ['at graph 0.96,1.17', 'maxrows 3', 'width -4', 'nobox', 'samplen 0.9'],
      lines = lines if noxerr else {},
      gpcalls = gpcalls,
      lmargin = 0.17, bmargin = 0.1, tmargin = 0.86, rmargin = 0.98,
      size = '9in,11in', arrow_offset = 0.9, #arrow_length = 0.4,
    )
  
    if nomed or noxerr or version == 'QM12': return 'done'
  
    # integrated enhancement factor
    if diffRel:
      enhance = {}
      data_enhance, medium_enhance, rhofo_enhance, vacrho_enhance = None, None, None, None
      for energy in sorted(data, key=float):
        for systLMR in [False, True]:
          suffix = str(energy)
          uEnhanceData = getMassRangesSums(
            data[energy], onlyLMR = True,
            systLMR = systLMR, suffix = suffix
          )
          uEnhanceCocktail = getMassRangesSums(
            cocktail[energy], onlyLMR = True,
            systLMR = systLMR, suffix = suffix
          )
          if energy in medium:
            uEnhanceMed = getMassRangesSums(
              medium[energy], onlyLMR = True,
              systLMR = systLMR, suffix = suffix
            )
          if energy in rhofo:
            uEnhanceRhoFO = getMassRangesSums(
              rhofo[energy], onlyLMR = True,
              systLMR = systLMR, suffix = suffix
            )
          if energy in vacrho:
            uEnhanceVacRho = getMassRangesSums(
              vacrho[energy], onlyLMR = True,
              systLMR = systLMR, suffix = suffix
            )
          if not systLMR: # uEnhance's are ufloats
            uEnhanceData /= uEnhanceCocktail
            dp = [
                float(energy), uEnhanceData.nominal_value, 0,
                getErrorComponent(uEnhanceData, 'stat'),
                getErrorComponent(uEnhanceData, 'syst')
            ]
            if data_enhance is None: data_enhance = [ dp ]
            else: data_enhance.append(dp)
            if energy in medium:
              uEnhanceMed /= uEnhanceCocktail
              dpM = [ float(energy), uEnhanceMed.nominal_value+1., 0, 0, 0 ]
              if medium_enhance is None: medium_enhance = [ dpM ]
              else: medium_enhance.append(dpM)
            if energy in rhofo:
              uEnhanceRhoFO /= uEnhanceCocktail
              dpM = [ float(energy), uEnhanceRhoFO.nominal_value+1., 0, 0, 0 ]
              if rhofo_enhance is None: rhofo_enhance = [ dpM ]
              else: rhofo_enhance.append(dpM)
            if energy in vacrho:
              uEnhanceVacRho /= uEnhanceCocktail
              dpM = [ float(energy), uEnhanceVacRho.nominal_value+1., 0, 0, 0 ]
              if vacrho_enhance is None: vacrho_enhance = [ dpM ]
              else: vacrho_enhance.append(dpM)
          else: # uEnhance's are dicts of ufloats
            for k in uEnhanceData:
              uEnhanceData[k] /= uEnhanceCocktail[k]
              dp = [
                float(energy), uEnhanceData[k].nominal_value, 0,
                getErrorComponent(uEnhanceData[k], 'stat'),
                getErrorComponent(uEnhanceData[k], 'syst')
              ]
              rngstr = k.split('_')[-1]
              data_key = 'data_' + rngstr
              if data_key not in enhance: enhance[data_key] = [ dp ]
              else: enhance[data_key].append(dp)
              if k in uEnhanceMed:
                uEnhanceMed[k] /= uEnhanceCocktail[k]
                dpM = [ float(energy), uEnhanceMed[k].nominal_value ]
                med_key = 'model_' + rngstr
                if med_key not in enhance: enhance[med_key] = [ dpM ]
                else: enhance[med_key].append(dpM)
              if k in uEnhanceRhoFO:
                uEnhanceRhoFO[k] /= uEnhanceCocktail[k]
                dpM = [ float(energy), uEnhanceRhoFO[k].nominal_value+1. ]
                rhofo_key = 'rhofo_' + rngstr
                if rhofo_key not in enhance: enhance[rhofo_key] = [ dpM ]
                else: enhance[rhofo_key].append(dpM)
              if k in uEnhanceVacRho:
                uEnhanceVacRho[k] /= uEnhanceCocktail[k]
                dpM = [ float(energy), uEnhanceVacRho[k].nominal_value+1. ]
                vacrho_key = 'vacrho_' + rngstr
                if vacrho_key not in enhance: enhance[vacrho_key] = [ dpM ]
                else: enhance[vacrho_key].append(dpM)
      xfacs = os.path.join(outDir, 'xfacs%s.dat' % version)
      if os.path.exists(xfacs): os.remove(xfacs)
      fSystLMR = open(xfacs, 'ab')
      for k in sorted(enhance.keys()):
        np.savetxt(fSystLMR, enhance[k], fmt = '%g', header = k, comments = '\n\n')
      fSystLMR.close()
      yr_upp = 4 if version == 'QM12Latest200' or version == 'QM14' else 7
      if version == 'LatestPatrickJieYi': yr_upp = 5.5
      #labels.update({
      #    '{LMR: %.2f < M_{ee} < %.2f GeV/c^{2}}' % (eRanges[1], eRanges[2]): [0.4,0.15,False]
      #})
      make_plot(
        data = [
            pseudo_point, pseudo_point, pseudo_point,
            np.array([[17.3,2.73,0,0.25,1.47]]),
            np.array([[200,4.7,0,0.4,1.5]]),
            np.array(enhance['data_0.15-0.75']),
            np.array(enhance['data_0.4-0.75']),
            np.array(medium_enhance),
            np.array(rhofo_enhance), np.array(vacrho_enhance)
        ],
        properties = [
            'lt 1 lw 4 ps 2 lc rgb "white" pt 19',
            'lt 1 lw 4 ps 2 lc rgb "white" pt 20',
            'lt 1 lw 4 ps 2 lc rgb "white" pt 18',
            'lt 1 lw 4 ps 2 lc %s pt 19' % default_colors[1],
            'lt 1 lw 4 ps 2 lc %s pt 20' % default_colors[3],
            'lt 1 lw 4 ps 2 lc %s pt 18' % default_colors[4],
            'lt 1 lw 4 ps 2 lc %s pt 18' % default_colors[0],
            'with lines lt 2 lw 4 lc %s' % default_colors[-1],
            'with lines lt 3 lw 4 lc %s' % default_colors[-1],
            'with lines lt 4 lw 4 lc %s' % default_colors[-1],
        ],
        titles = [
            'CERES Pb+Au', 'PHENIX Au+Au', 'STAR Au+Au',
            '', '', '', '',
            'HMBT + QGP', 'BW/FO-{/Symbol \162}', '{/Symbol \162}/{/Symbol \167} VacSF+FB',
        ],
        name = os.path.join(outDir, 'enhance%s' % version),
        xlabel = '{/Symbol \326}s_{NN} (GeV)',
        ylabel = 'LMR Enhancement Factor',
        xlog = True, key = [ 'at graph 0.9,0.98', 'nobox', 'maxrows 4' ],
        size = '10in,8in', bmargin = 0.13, tmargin = 0.92, rmargin = 0.99,
        yr = [1.,yr_upp], xr = [14,220], gpcalls = [
          'format x "%g"',
          'xtics (10, 20,"" 30, 40,"" 50, 60,"" 70,"" 80,"" 90, 100, 200)',
          'boxwidth 0.025 absolute',
          'label 50 "{/=18 0.2 < M_{ee} < 0.6 GeV/c^{2}}" at 15.5,3 tc %s rotate center' % default_colors[1],
          'label 51 "{/=18 0.15 < M_{ee} < 0.75 GeV/c^{2}}" at 180,4.2 tc %s rotate center' % default_colors[3],
          'label 52 "{/=18 0.4 < M_{ee} < 0.75 GeV/c^{2}}" at 75,3.1 tc %s rotate by -20' % default_colors[0],
          'label 53 "{/=18 0.15 < M_{ee} < 0.75 GeV/c^{2}}" at 50,1.2 tc %s' % default_colors[4]
        ], #labels = labels
      )
      return 'done'
  
    # integrated excess yield in mass ranges
    excess = {}
    for k, v in dataOrdered.iteritems():
      suffix = ''
      energy = getEnergy4Key(re.compile('\d+').search(k).group())
      if fnmatch(k, '*Med.*'):
        suffix = '_Med'
        if version != 'LatestPatrickJieYi' and energy == '27': continue # TODO
      if fnmatch(k, '*RhoFO.*'): suffix = '_RhoFO'
      if fnmatch(k, '*VacRho.*'): suffix = '_VacRho'
      exc = getMassRangesSums(np.array(v), onlyLMR = True)
      if divdNdy: exc /= dNdyPi0[energy] * 1e-2
      dp = [
          float(energy), exc.nominal_value, 0,
          getErrorComponent(exc, 'stat'), getErrorComponent(exc, 'syst')
      ]
      if suffix == '_Med' and not diffRel and not divdNdy:
          print dp
      key = 'LMR' + suffix
      if key not in excess: excess[key] = [ dp ]
      else: excess[key].append(dp)
    logging.debug(excess)
    avdata = np.array(excess['LMR'])
    avg = np.average(avdata[:,1], weights = avdata[:,3])
    graph_data = [
        np.array([
            [ 7.7, avg, 0, 0, avdata[-1][-1]],
            [ 19.6, avg, 0, 0, avdata[-1][-1]]
        ]),
        np.array([
            [ 19.6, avg, 0, 0, 0], [ 200., avg, 0, 0, 0]
        ]),
        np.array([
            [ 7.7, 2*avg, 0, 0, 0], [ 19.6, avg, 0, 0, 0],
        ]),
        np.array(excess['LMR']),
    ]
    props = [
        'with filledcurves pt 0 lc %s lw 4 lt 2' % default_colors[8],
        'with lines lc %s lw 4 lt 2' % default_colors[8],
        'with lines lc %s lw 8 lt 2' % default_colors[1],
        'lt 1 lw 4 ps 2 lc %s pt 18' % default_colors[0],
    ]
    tits = [
        'BES-I extrapolation', '', 'model expectation', 'STAR Au+Au',
    ]
    if version != 'QM14':
        graph_data += [
            np.array(excess['LMR_Med']),
            np.array(excess['LMR_VacRho']),
            np.array(excess['LMR_RhoFO']),
        ]
        props += [
            'with lines lt 2 lw 4 lc %s' % default_colors[-1],
            'with lines lt 3 lw 4 lc %s' % default_colors[-1],
            'with lines lt 4 lw 4 lc %s' % default_colors[-1],
        ]
        tits += [
            'HMBT + QGP', '{/Symbol \162}/{/Symbol \167} VacSF+FB', 'BW/FO-{/Symbol \162}',
        ]
    yr_upp = 4.5 if version == 'QM12Latest200' or version == 'QM14' else 7
    if version == 'LatestPatrickJieYi': yr_upp = 2 if divdNdy else 2.
    labels = {} if version != 'QM14' else labels
    if divdNdy:
      labels.update(dict((str(v), [float(k)*0.9,yr_upp*1.05,True]) for k,v in dNdyPi0.items()))
      labels.update({ 'dN/dy|_{/Symbol \\160}': [100,yr_upp*1.05,True]})
    gpcalls = [
      'format x "%g"',
      'xtics (7,10,20,"" 30, 40,"" 50, 60,"" 70,"" 80,"" 90, 100, 200)',
      'boxwidth 0.025 absolute',
    ]
    if version == 'QM14':
      labels.update({
        '{LMR: %.2f < M_{ee} < %.2f GeV/c^{2}}' % (eRanges[1], eRanges[2]): [0.4,0.15,False],
      })
    else:
      gpcalls.append('label 52 "{/=18 0.4 < M_{ee} < 0.75 GeV/c^{2}}" at 60,0.4 tc %s' % default_colors[0])
    make_plot(
      data = graph_data, properties = props, titles = tits,
      name = os.path.join(outDir, 'excess%s%s' % (version,'DivdNdy' if divdNdy else '')),
      xlabel = '{/Symbol \326}s_{NN} (GeV)',
      ylabel = 'LMR Excess Yield %s({/Symbol \264} 10^{-%d})' % (
          '/ dN/dy|_{/Symbol \\160}  ' if divdNdy else '', 5 if divdNdy else 3
      ),
      xlog = True, xr = [7,220], size = '10in,8in',
      key = ['at graph 1.05,0.98', 'width -3', 'nobox', 'maxrows 3'],
      bmargin = 0.13, tmargin = 0.92, rmargin = 0.99,
      yr = [0,yr_upp], gpcalls = gpcalls, labels = labels,
    )
    return 'done'