gp_ccX: first version done

pyana/examples/gp_ccX.py

@@ -5,6 +5,7 @@
 from ..ccsgp.config import default_colors
 from scipy.optimize import curve_fit
 import numpy as np
+from uncertainties import ufloat
 
 def linear(x, a, b):
     return a*x+b
@@ -43,22 +44,49 @@ def gp_ccX():
         lambda x, c, d: fitfunc(x, popt1[0], popt1[1], c, d),
         alldata[:,0], alldata[:,1], sigma=alldata[:,3], absolute_sigma=True,
     )
-    fitdata = np.array([[
-        x, fitfunc(x, popt1[0], popt1[1], popt2[0], popt2[1]), 0, 0, 0
-    ] for x in np.linspace(1,4)])
+    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 = np.sqrt(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 = data.values() + [fitdata],
+        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 2'],
-        titles = data.keys() + ['parameterization'],
-        xlabel = 'log_{10}[{/Symbol \326}s_{NN} (GeV)]',
-        ylabel = 'log_{10}[{/Symbol \163}@_{c@^{/=18-}c}^{NN} ({/Symbol \155}b)]',
+        ] + ['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 = '8.5in,8in', xr = [1, 4], yr = [0.5,4.3],
-        key = ['bottom right', 'nobox', 'width -3'],
-        tmargin = 0.99, rmargin = 0.97, bmargin = 0.13
+        size = '8.5in,8in', 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,
     )
 
 if __name__ == '__main__':