gp_rdiff: calc excess yields + errors

tschaume · Mar 13, 2014 · faf69c9 · faf69c9
1 parent d53e55b
commit faf69c9
Show file tree

Hide file tree

Showing 2 changed files with 22 additions and 9 deletions.
diff --git a/pyana/examples/gp_rdiff.py b/pyana/examples/gp_rdiff.py
@@ -3,7 +3,7 @@
 from fnmatch import fnmatch
 from collections import OrderedDict
 from .utils import getWorkDirs, checkSymLink
-from .utils import getUArray, getEdges, getCocktailSum, enumzipEdges
+from .utils import getUArray, getEdges, getCocktailSum, enumzipEdges, getMassRangesSums
 from ..ccsgp.ccsgp import make_plot
 from ..ccsgp.utils import getOpts, zip_flat
 from ..ccsgp.config import default_colors
@@ -61,13 +61,14 @@ def gp_rdiff(version):
         uCocktailSum = getCocktailSum(e0, e1, eCocktail, uCocktail)
         # calc. difference and divide by data binwidth again
         # + set data point
-        xs = xshift if energy == '39' else 0.
+        #xs = xshift if energy == '39' else 0.
+        xs = 0. # TODO: cannot use shift when calculating excess yield based on dataOrdered
         if not l:
           uDiff = uData[i] - uCocktailSum
           uDiff /= data[energy][i,2] * 2 * yunit
           dp = [
             data[energy][i,0] + xs, uDiff.nominal_value,
-            0., data[energy][i,3] / yunit, uDiff.std_dev
+            data[energy][i,2], data[energy][i,3] / yunit, uDiff.std_dev
           ]
           key = ' '.join([energy, 'GeV'])
         else:
@@ -85,6 +86,14 @@ def gp_rdiff(version):
         if key in dataOrdered: dataOrdered[key].append(dp)
         else: dataOrdered[key] = [ dp ]
 
+  # integrated excess yield in mass ranges
+  excess = {}
+  for k, v in dataOrdered.iteritems():
+    if fnmatch(k, '*Med.*'): continue
+    energy = re.compile('\d+').search(k).group()
+    getMassRangesSums(energy, np.array(v), excess, onlyLMR = True)
+  print excess
+
   # make plot
   nSets = len(dataOrdered)
   nSetsPlot = nSets/2 if nSets > 4 else nSets

diff --git a/pyana/examples/utils.py b/pyana/examples/utils.py
@@ -35,6 +35,7 @@ def getWorkDirs():
 
 def getUArray(npArr):
   """uncertainty array multiplied by binwidth (col2 = dx)"""
+  # propagates systematic uncertainty!
   return unp.uarray(npArr[:,1], npArr[:,4]) * npArr[:,2] * 2
 
 def getEdges(npArr):
@@ -98,14 +99,17 @@ def getCocktailSum(e0, e1, eCocktail, uCocktail):
     logging.debug('    sum: {}'.format(uCocktailSum))
   return uCocktailSum
 
-def getMassRangesSums(energy, indata, outdata):
-  uMedium = getUArray(indata)
-  eMedium = getEdges(indata)
+def getMassRangesSums(energy, indata, outdata, onlyLMR = False):
+  # combine stat. and syst. errorbars
+  indata[:,4] = np.sqrt(indata[:,3] * indata[:,3] + indata[:,4] * indata[:,4])
+  uInData = getUArray(indata)
+  eInData = getEdges(indata)
   for i, (e0, e1) in enumzipEdges(eRanges):
-    uSum = getCocktailSum(e0, e1, eMedium, uMedium)
+    if onlyLMR and i != 1: continue
+    uSum = getCocktailSum(e0, e1, eInData, uInData)
     logging.debug('%s> %g - %g: %r' % (energy, e0, e1, uSum))
     dp = [
-      float(energy), uSum.nominal_value, 0, uSum.std_dev, 0
-    ] # TODO: syst.  uncertainties
+      float(energy), uSum.nominal_value, 0, 0, uSum.std_dev
+    ]
     if mass_titles[i] in outdata: outdata[mass_titles[i]].append(dp)
     else: outdata[mass_titles[i]] = [ dp ]