From d85bff0d81ac03f4eab7efd86dc5ff0e4c78f5c9 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Tue, 24 Sep 2013 23:23:50 -0700
Subject: [PATCH 01/23] initial draft (spectrum similar to fbu preprint)

---
 python/generate_bump.py | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)
 create mode 100755 python/generate_bump.py

diff --git a/python/generate_bump.py b/python/generate_bump.py
new file mode 100755
index 0000000..3d57611
--- /dev/null
+++ b/python/generate_bump.py
@@ -0,0 +1,23 @@
+#!/bin/env python
+
+import numpy as np
+import matplotlib.mlab as mlab
+import matplotlib.pyplot as plt
+from math import sqrt
+
+offset, scale = 100.0, 500.0
+num_bins = 50
+x = np.random.exponential(scale=1.0, size=1e6) * scale + offset
+xm = np.random.normal(loc=4.0, scale=0.25, size=1e5) * scale + offset
+xtot = np.concatenate((x, xm))
+
+def smear(x, a=0.5, b=0.1) :
+   sigma = a*sqrt(x) + b*x
+   return x+np.random.normal(loc=0.0, scale=sigma)
+
+xsmear = [smear(x) for x in xtot]
+n, bins, patches = plt.hist(xtot, num_bins, facecolor='green', alpha=0.25, log=True)
+n, bins, patches = plt.hist(xsmear, num_bins, facecolor='blue', alpha=0.25, log=True)
+
+plt.show()
+plt.savefig('fallingExp_withBump.png')

From 48c010216e1cc879b128bc83f8f053ed61030009 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Mon, 7 Oct 2013 00:59:25 -0700
Subject: [PATCH 02/23] significant re-write

Details:
- factorize the code in well-defined functions
- add docstrings explaining what I am doing
- rewrite most of this script. Now plotting and using the resulting
histograms as bins/counts that can be used for the unfolding.
- compute normalized response matrix
- add plotting functions
---
 python/generate_bump.py | 61 +++++++++++++++++++++++++++++++++--------
 1 file changed, 49 insertions(+), 12 deletions(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 3d57611..892044c 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -5,19 +5,56 @@
 import matplotlib.pyplot as plt
 from math import sqrt
 
-offset, scale = 100.0, 500.0
-num_bins = 50
-x = np.random.exponential(scale=1.0, size=1e6) * scale + offset
-xm = np.random.normal(loc=4.0, scale=0.25, size=1e5) * scale + offset
-xtot = np.concatenate((x, xm))
+def generateTruthSpectrum(offset=100.0, scale=500.0, nbkg=1e6, nsig=1e5) :
+    """Generate pseudo-events with an observable in [offset, scale].
+    The observable spectrum is inspired to the one expected for a
+    di-jet resonance, where the background is distributed as a falling
+    exponential, and the signal is a gaussian on top of it.
+    """
+    bkg = np.random.exponential(scale=1.0, size=nbkg) * scale + offset
+    sig = np.random.normal(loc=4.0, scale=0.25, size=nsig) * scale + offset
+    tot = np.concatenate((bkg, sig))
+    return tot
 
 def smear(x, a=0.5, b=0.1) :
-   sigma = a*sqrt(x) + b*x
-   return x+np.random.normal(loc=0.0, scale=sigma)
+    """Smear the pseudo-events generated with generateTruthSpectrum as
+    to mimic the jet resolution, i.e. with a gaussianely-distributed
+    \delta, see eq. 13, eq. 14, and par. 7.1 of the FBU paper.
+    """
+    sigma = a*sqrt(x) + b*x
+    return x+np.random.normal(loc=0.0, scale=sigma)
 
-xsmear = [smear(x) for x in xtot]
-n, bins, patches = plt.hist(xtot, num_bins, facecolor='green', alpha=0.25, log=True)
-n, bins, patches = plt.hist(xsmear, num_bins, facecolor='blue', alpha=0.25, log=True)
+sigPlusBkg = generateTruthSpectrum()
+spbSmeared = [smear(x) for x in sigPlusBkg]
+numBins = 50
 
-plt.show()
-plt.savefig('fallingExp_withBump.png')
+nBinsTruth, truthMin, truthMax = 50, 500.0, 4500.0
+nBinsReco,  recoMin,  recoMax  = 50, 500.0, 4500.0
+histTruth, binsTruth        = np.histogram(sigPlusBkg, bins=nBinsTruth, range=(truthMin, truthMax))
+histReco,  binsReco         = np.histogram(spbSmeared, bins=nBinsReco, range=(recoMin, recoMax))
+respHist,  xedges, yedges   = np.histogram2d(sigPlusBkg, spbSmeared, bins=(nBinsTruth, nBinsReco),
+                                             range=((truthMin, truthMax), (recoMin,  recoMax)))
+def normalized(mat) :
+    mat = mat.astype(float)
+    return mat / mat.sum()
+respMat = normalized(respHist)
+print respHist
+print respMat
+
+def plotTruthAndReco(histT, binsT, histR, binsR, outfname) :
+    def getWidth(bins) : return bins[1] - bins[0]
+    plt.figure()
+    plt.bar(binsT[:-1], histT, color='b', width=getWidth(binsT), alpha=0.15)
+    plt.bar(binsR[:-1], histR, color='r', width=getWidth(binsR), alpha=0.15)
+    plt.savefig(outfname)
+
+def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
+    from matplotlib.colors import LogNorm
+    extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
+    plt.figure()
+    plt.imshow(respMatrix, extent=extent, interpolation='nearest', norm=LogNorm(), origin='lower')
+    plt.colorbar()
+    plt.savefig(outfname)
+
+plotTruthAndReco(histTruth, binsTruth, histReco, binsReco, 'fallingExp_withBump.png')
+plotRespMatrix(respHist, xedges, yedges, 'responseMatrix.png')

From 49722d798019bc3a2546c91193876ae0813114ab Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Mon, 7 Oct 2013 01:21:10 -0700
Subject: [PATCH 03/23] utils to write/read np.array from/to json

---
 python/utils.py | 10 ++++++++++
 1 file changed, 10 insertions(+)
 create mode 100644 python/utils.py

diff --git a/python/utils.py b/python/utils.py
new file mode 100644
index 0000000..e409926
--- /dev/null
+++ b/python/utils.py
@@ -0,0 +1,10 @@
+import json
+import numpy as np
+
+def array2json(a, outfname) :
+    with open(outfname, 'w') as out:
+        json.dump(a.tolist(), out)
+
+def json2array(infname) :
+    with open(infname) as inp:
+        return numpy.array(json.load(inp))

From 72d42d271d17a22277b59da47b0aacd05a0225a5 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Mon, 7 Oct 2013 01:42:11 -0700
Subject: [PATCH 04/23] separate the generation from the unfolding

Details:
two separate logic steps: generate the pseudo-data, and then unfold.
The pseudo-data can come either from the json file, or be generated.
---
 python/generate_bump.py | 80 +++++++++++++++++++++++++----------------
 1 file changed, 50 insertions(+), 30 deletions(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 892044c..2b4c6b9 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -4,6 +4,11 @@
 import matplotlib.mlab as mlab
 import matplotlib.pyplot as plt
 from math import sqrt
+import os
+from utils import array2json, json2array
+
+
+regenerate = True
 
 def generateTruthSpectrum(offset=100.0, scale=500.0, nbkg=1e6, nsig=1e5) :
     """Generate pseudo-events with an observable in [offset, scale].
@@ -24,37 +29,52 @@ def smear(x, a=0.5, b=0.1) :
     sigma = a*sqrt(x) + b*x
     return x+np.random.normal(loc=0.0, scale=sigma)
 
-sigPlusBkg = generateTruthSpectrum()
-spbSmeared = [smear(x) for x in sigPlusBkg]
-numBins = 50
+jsonDir = 'data/bump/'
+jsonTruthFname  = jsonDir+'truth.json'
+jsonSmearFname  = jsonDir+'reco.json'
+jsonResMatFname = jsonDir+'resMat.json'
 
-nBinsTruth, truthMin, truthMax = 50, 500.0, 4500.0
-nBinsReco,  recoMin,  recoMax  = 50, 500.0, 4500.0
-histTruth, binsTruth        = np.histogram(sigPlusBkg, bins=nBinsTruth, range=(truthMin, truthMax))
-histReco,  binsReco         = np.histogram(spbSmeared, bins=nBinsReco, range=(recoMin, recoMax))
-respHist,  xedges, yedges   = np.histogram2d(sigPlusBkg, spbSmeared, bins=(nBinsTruth, nBinsReco),
-                                             range=((truthMin, truthMax), (recoMin,  recoMax)))
 def normalized(mat) :
     mat = mat.astype(float)
     return mat / mat.sum()
-respMat = normalized(respHist)
-print respHist
-print respMat
-
-def plotTruthAndReco(histT, binsT, histR, binsR, outfname) :
-    def getWidth(bins) : return bins[1] - bins[0]
-    plt.figure()
-    plt.bar(binsT[:-1], histT, color='b', width=getWidth(binsT), alpha=0.15)
-    plt.bar(binsR[:-1], histR, color='r', width=getWidth(binsR), alpha=0.15)
-    plt.savefig(outfname)
-
-def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
-    from matplotlib.colors import LogNorm
-    extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
-    plt.figure()
-    plt.imshow(respMatrix, extent=extent, interpolation='nearest', norm=LogNorm(), origin='lower')
-    plt.colorbar()
-    plt.savefig(outfname)
-
-plotTruthAndReco(histTruth, binsTruth, histReco, binsReco, 'fallingExp_withBump.png')
-plotRespMatrix(respHist, xedges, yedges, 'responseMatrix.png')
+
+def generateAndPlot() :
+    sigPlusBkg = generateTruthSpectrum()
+    spbSmeared = np.array([smear(x) for x in sigPlusBkg])
+    numBins = 50
+    nBinsTruth, truthMin, truthMax = 50, 500.0, 4500.0
+    nBinsReco,  recoMin,  recoMax  = 50, 500.0, 4500.0
+    histTruth, binsTruth        = np.histogram(sigPlusBkg, bins=nBinsTruth, range=(truthMin, truthMax))
+    histReco,  binsReco         = np.histogram(spbSmeared, bins=nBinsReco, range=(recoMin, recoMax))
+    respHist,  xedges, yedges   = np.histogram2d(sigPlusBkg, spbSmeared, bins=(nBinsTruth, nBinsReco),
+                                                 range=((truthMin, truthMax), (recoMin,  recoMax)))
+    respMat = normalized(respHist)
+    def plotTruthAndReco(histT, binsT, histR, binsR, outfname) :
+        def getWidth(bins) : return bins[1] - bins[0]
+        plt.figure()
+        plt.bar(binsT[:-1], histT, color='b', width=getWidth(binsT), alpha=0.15)
+        plt.bar(binsR[:-1], histR, color='r', width=getWidth(binsR), alpha=0.15)
+        plt.savefig(outfname)
+    def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
+        from matplotlib.colors import LogNorm
+        extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
+        plt.figure()
+        plt.imshow(respMatrix, extent=extent, interpolation='nearest', norm=LogNorm(), origin='lower')
+        plt.colorbar()
+        plt.savefig(outfname)
+    plotTruthAndReco(histTruth, binsTruth, histReco, binsReco, 'fallingExp_withBump.png')
+    plotRespMatrix(respHist, xedges, yedges, 'responseMatrix.png')
+    return histTruth, histReco, respMat
+
+histTruth, histReco, respMat = None, None, None
+if regenerate :
+    histTruth, histReco, respMat = generateAndPlot()
+    if not os.path.isdir(jsonDir) : os.path.mkdir(jsonDir)
+    for a,f in [(histTruth, jsonTruthFname),
+                (histReco,  jsonSmearFname),
+                (respMat,   jsonResMatFname)] :
+        array2json(a,f)
+else :
+    histTruth  = json2array(jsonTruthFname)
+    histReco   = json2array(jsonSmearFname),
+    respMat    = json2array(jsonResMatFname)

From e59f6321c9f2a5a6d0d53909de8d74ae3fc91ae2 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Mon, 7 Oct 2013 01:46:46 -0700
Subject: [PATCH 05/23] almost there, but pyfbu does not work

Details:
complaining about template. I thought it wasn't needed anymore.
Need to ask clement and Francesco about it.
---
 python/generate_bump.py | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 2b4c6b9..18d0b47 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -6,7 +6,7 @@
 from math import sqrt
 import os
 from utils import array2json, json2array
-
+from pyFBU import pyFBU
 
 regenerate = True
 
@@ -78,3 +78,20 @@ def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
     histTruth  = json2array(jsonTruthFname)
     histReco   = json2array(jsonSmearFname),
     respMat    = json2array(jsonResMatFname)
+
+pyfbu = pyFBU()
+pyfbu.nMCMC        = 100000
+pyfbu.nBurn        = 1000
+pyfbu.nThin        = 10
+pyfbu.lower        = 70000
+pyfbu.upper        = 140000
+pyfbu.jsonData     = jsonSmearFname
+pyfbu.jsonMig      = jsonResMatFname
+pyfbu.jsonBkg      = '' # this shouldnt be a requirement; pass [0.] ?
+pyfbu.templateFile = '' # templateFile  ?
+pyfbu.modelName    = '' # modelName if modelName else pyfbu.modelName ? I thought this wasn't necessary anymore
+pyfbu.verbose      = True # verbose ? will be on cmd-line
+
+pyfbu.run()
+
+trace = pyfbu.trace

From 99ecb5c97d9cff664aa1b44ba57e533d3bed3240 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Tue, 24 Sep 2013 23:23:50 -0700
Subject: [PATCH 06/23] initial draft (spectrum similar to fbu preprint)

---
 python/generate_bump.py | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)
 create mode 100755 python/generate_bump.py

diff --git a/python/generate_bump.py b/python/generate_bump.py
new file mode 100755
index 0000000..3d57611
--- /dev/null
+++ b/python/generate_bump.py
@@ -0,0 +1,23 @@
+#!/bin/env python
+
+import numpy as np
+import matplotlib.mlab as mlab
+import matplotlib.pyplot as plt
+from math import sqrt
+
+offset, scale = 100.0, 500.0
+num_bins = 50
+x = np.random.exponential(scale=1.0, size=1e6) * scale + offset
+xm = np.random.normal(loc=4.0, scale=0.25, size=1e5) * scale + offset
+xtot = np.concatenate((x, xm))
+
+def smear(x, a=0.5, b=0.1) :
+   sigma = a*sqrt(x) + b*x
+   return x+np.random.normal(loc=0.0, scale=sigma)
+
+xsmear = [smear(x) for x in xtot]
+n, bins, patches = plt.hist(xtot, num_bins, facecolor='green', alpha=0.25, log=True)
+n, bins, patches = plt.hist(xsmear, num_bins, facecolor='blue', alpha=0.25, log=True)
+
+plt.show()
+plt.savefig('fallingExp_withBump.png')

From 9eb363838e0fd2ead4d996c2eb0b7bbb8a55a6b6 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Mon, 7 Oct 2013 00:59:25 -0700
Subject: [PATCH 07/23] significant re-write

Details:
- factorize the code in well-defined functions
- add docstrings explaining what I am doing
- rewrite most of this script. Now plotting and using the resulting
histograms as bins/counts that can be used for the unfolding.
- compute normalized response matrix
- add plotting functions
---
 python/generate_bump.py | 61 +++++++++++++++++++++++++++++++++--------
 1 file changed, 49 insertions(+), 12 deletions(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 3d57611..892044c 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -5,19 +5,56 @@
 import matplotlib.pyplot as plt
 from math import sqrt
 
-offset, scale = 100.0, 500.0
-num_bins = 50
-x = np.random.exponential(scale=1.0, size=1e6) * scale + offset
-xm = np.random.normal(loc=4.0, scale=0.25, size=1e5) * scale + offset
-xtot = np.concatenate((x, xm))
+def generateTruthSpectrum(offset=100.0, scale=500.0, nbkg=1e6, nsig=1e5) :
+    """Generate pseudo-events with an observable in [offset, scale].
+    The observable spectrum is inspired to the one expected for a
+    di-jet resonance, where the background is distributed as a falling
+    exponential, and the signal is a gaussian on top of it.
+    """
+    bkg = np.random.exponential(scale=1.0, size=nbkg) * scale + offset
+    sig = np.random.normal(loc=4.0, scale=0.25, size=nsig) * scale + offset
+    tot = np.concatenate((bkg, sig))
+    return tot
 
 def smear(x, a=0.5, b=0.1) :
-   sigma = a*sqrt(x) + b*x
-   return x+np.random.normal(loc=0.0, scale=sigma)
+    """Smear the pseudo-events generated with generateTruthSpectrum as
+    to mimic the jet resolution, i.e. with a gaussianely-distributed
+    \delta, see eq. 13, eq. 14, and par. 7.1 of the FBU paper.
+    """
+    sigma = a*sqrt(x) + b*x
+    return x+np.random.normal(loc=0.0, scale=sigma)
 
-xsmear = [smear(x) for x in xtot]
-n, bins, patches = plt.hist(xtot, num_bins, facecolor='green', alpha=0.25, log=True)
-n, bins, patches = plt.hist(xsmear, num_bins, facecolor='blue', alpha=0.25, log=True)
+sigPlusBkg = generateTruthSpectrum()
+spbSmeared = [smear(x) for x in sigPlusBkg]
+numBins = 50
 
-plt.show()
-plt.savefig('fallingExp_withBump.png')
+nBinsTruth, truthMin, truthMax = 50, 500.0, 4500.0
+nBinsReco,  recoMin,  recoMax  = 50, 500.0, 4500.0
+histTruth, binsTruth        = np.histogram(sigPlusBkg, bins=nBinsTruth, range=(truthMin, truthMax))
+histReco,  binsReco         = np.histogram(spbSmeared, bins=nBinsReco, range=(recoMin, recoMax))
+respHist,  xedges, yedges   = np.histogram2d(sigPlusBkg, spbSmeared, bins=(nBinsTruth, nBinsReco),
+                                             range=((truthMin, truthMax), (recoMin,  recoMax)))
+def normalized(mat) :
+    mat = mat.astype(float)
+    return mat / mat.sum()
+respMat = normalized(respHist)
+print respHist
+print respMat
+
+def plotTruthAndReco(histT, binsT, histR, binsR, outfname) :
+    def getWidth(bins) : return bins[1] - bins[0]
+    plt.figure()
+    plt.bar(binsT[:-1], histT, color='b', width=getWidth(binsT), alpha=0.15)
+    plt.bar(binsR[:-1], histR, color='r', width=getWidth(binsR), alpha=0.15)
+    plt.savefig(outfname)
+
+def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
+    from matplotlib.colors import LogNorm
+    extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
+    plt.figure()
+    plt.imshow(respMatrix, extent=extent, interpolation='nearest', norm=LogNorm(), origin='lower')
+    plt.colorbar()
+    plt.savefig(outfname)
+
+plotTruthAndReco(histTruth, binsTruth, histReco, binsReco, 'fallingExp_withBump.png')
+plotRespMatrix(respHist, xedges, yedges, 'responseMatrix.png')

From c2a1f51385e8075a00f660e111f7114008444862 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Mon, 7 Oct 2013 01:21:10 -0700
Subject: [PATCH 08/23] utils to write/read np.array from/to json

---
 python/utils.py | 10 ++++++++++
 1 file changed, 10 insertions(+)
 create mode 100644 python/utils.py

diff --git a/python/utils.py b/python/utils.py
new file mode 100644
index 0000000..e409926
--- /dev/null
+++ b/python/utils.py
@@ -0,0 +1,10 @@
+import json
+import numpy as np
+
+def array2json(a, outfname) :
+    with open(outfname, 'w') as out:
+        json.dump(a.tolist(), out)
+
+def json2array(infname) :
+    with open(infname) as inp:
+        return numpy.array(json.load(inp))

From 02fa5f47ab070acd4f823b85cd9ffc22ce86d73a Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Mon, 7 Oct 2013 01:42:11 -0700
Subject: [PATCH 09/23] separate the generation from the unfolding

Details:
two separate logic steps: generate the pseudo-data, and then unfold.
The pseudo-data can come either from the json file, or be generated.
---
 python/generate_bump.py | 80 +++++++++++++++++++++++++----------------
 1 file changed, 50 insertions(+), 30 deletions(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 892044c..2b4c6b9 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -4,6 +4,11 @@
 import matplotlib.mlab as mlab
 import matplotlib.pyplot as plt
 from math import sqrt
+import os
+from utils import array2json, json2array
+
+
+regenerate = True
 
 def generateTruthSpectrum(offset=100.0, scale=500.0, nbkg=1e6, nsig=1e5) :
     """Generate pseudo-events with an observable in [offset, scale].
@@ -24,37 +29,52 @@ def smear(x, a=0.5, b=0.1) :
     sigma = a*sqrt(x) + b*x
     return x+np.random.normal(loc=0.0, scale=sigma)
 
-sigPlusBkg = generateTruthSpectrum()
-spbSmeared = [smear(x) for x in sigPlusBkg]
-numBins = 50
+jsonDir = 'data/bump/'
+jsonTruthFname  = jsonDir+'truth.json'
+jsonSmearFname  = jsonDir+'reco.json'
+jsonResMatFname = jsonDir+'resMat.json'
 
-nBinsTruth, truthMin, truthMax = 50, 500.0, 4500.0
-nBinsReco,  recoMin,  recoMax  = 50, 500.0, 4500.0
-histTruth, binsTruth        = np.histogram(sigPlusBkg, bins=nBinsTruth, range=(truthMin, truthMax))
-histReco,  binsReco         = np.histogram(spbSmeared, bins=nBinsReco, range=(recoMin, recoMax))
-respHist,  xedges, yedges   = np.histogram2d(sigPlusBkg, spbSmeared, bins=(nBinsTruth, nBinsReco),
-                                             range=((truthMin, truthMax), (recoMin,  recoMax)))
 def normalized(mat) :
     mat = mat.astype(float)
     return mat / mat.sum()
-respMat = normalized(respHist)
-print respHist
-print respMat
-
-def plotTruthAndReco(histT, binsT, histR, binsR, outfname) :
-    def getWidth(bins) : return bins[1] - bins[0]
-    plt.figure()
-    plt.bar(binsT[:-1], histT, color='b', width=getWidth(binsT), alpha=0.15)
-    plt.bar(binsR[:-1], histR, color='r', width=getWidth(binsR), alpha=0.15)
-    plt.savefig(outfname)
-
-def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
-    from matplotlib.colors import LogNorm
-    extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
-    plt.figure()
-    plt.imshow(respMatrix, extent=extent, interpolation='nearest', norm=LogNorm(), origin='lower')
-    plt.colorbar()
-    plt.savefig(outfname)
-
-plotTruthAndReco(histTruth, binsTruth, histReco, binsReco, 'fallingExp_withBump.png')
-plotRespMatrix(respHist, xedges, yedges, 'responseMatrix.png')
+
+def generateAndPlot() :
+    sigPlusBkg = generateTruthSpectrum()
+    spbSmeared = np.array([smear(x) for x in sigPlusBkg])
+    numBins = 50
+    nBinsTruth, truthMin, truthMax = 50, 500.0, 4500.0
+    nBinsReco,  recoMin,  recoMax  = 50, 500.0, 4500.0
+    histTruth, binsTruth        = np.histogram(sigPlusBkg, bins=nBinsTruth, range=(truthMin, truthMax))
+    histReco,  binsReco         = np.histogram(spbSmeared, bins=nBinsReco, range=(recoMin, recoMax))
+    respHist,  xedges, yedges   = np.histogram2d(sigPlusBkg, spbSmeared, bins=(nBinsTruth, nBinsReco),
+                                                 range=((truthMin, truthMax), (recoMin,  recoMax)))
+    respMat = normalized(respHist)
+    def plotTruthAndReco(histT, binsT, histR, binsR, outfname) :
+        def getWidth(bins) : return bins[1] - bins[0]
+        plt.figure()
+        plt.bar(binsT[:-1], histT, color='b', width=getWidth(binsT), alpha=0.15)
+        plt.bar(binsR[:-1], histR, color='r', width=getWidth(binsR), alpha=0.15)
+        plt.savefig(outfname)
+    def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
+        from matplotlib.colors import LogNorm
+        extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
+        plt.figure()
+        plt.imshow(respMatrix, extent=extent, interpolation='nearest', norm=LogNorm(), origin='lower')
+        plt.colorbar()
+        plt.savefig(outfname)
+    plotTruthAndReco(histTruth, binsTruth, histReco, binsReco, 'fallingExp_withBump.png')
+    plotRespMatrix(respHist, xedges, yedges, 'responseMatrix.png')
+    return histTruth, histReco, respMat
+
+histTruth, histReco, respMat = None, None, None
+if regenerate :
+    histTruth, histReco, respMat = generateAndPlot()
+    if not os.path.isdir(jsonDir) : os.path.mkdir(jsonDir)
+    for a,f in [(histTruth, jsonTruthFname),
+                (histReco,  jsonSmearFname),
+                (respMat,   jsonResMatFname)] :
+        array2json(a,f)
+else :
+    histTruth  = json2array(jsonTruthFname)
+    histReco   = json2array(jsonSmearFname),
+    respMat    = json2array(jsonResMatFname)

From ba2a605f21fa0d15d677b3141ab98a8519192736 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Mon, 7 Oct 2013 01:46:46 -0700
Subject: [PATCH 10/23] almost there, but pyfbu does not work

Details:
complaining about template. I thought it wasn't needed anymore.
Need to ask clement and Francesco about it.
---
 python/generate_bump.py | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 2b4c6b9..18d0b47 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -6,7 +6,7 @@
 from math import sqrt
 import os
 from utils import array2json, json2array
-
+from pyFBU import pyFBU
 
 regenerate = True
 
@@ -78,3 +78,20 @@ def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
     histTruth  = json2array(jsonTruthFname)
     histReco   = json2array(jsonSmearFname),
     respMat    = json2array(jsonResMatFname)
+
+pyfbu = pyFBU()
+pyfbu.nMCMC        = 100000
+pyfbu.nBurn        = 1000
+pyfbu.nThin        = 10
+pyfbu.lower        = 70000
+pyfbu.upper        = 140000
+pyfbu.jsonData     = jsonSmearFname
+pyfbu.jsonMig      = jsonResMatFname
+pyfbu.jsonBkg      = '' # this shouldnt be a requirement; pass [0.] ?
+pyfbu.templateFile = '' # templateFile  ?
+pyfbu.modelName    = '' # modelName if modelName else pyfbu.modelName ? I thought this wasn't necessary anymore
+pyfbu.verbose      = True # verbose ? will be on cmd-line
+
+pyfbu.run()
+
+trace = pyfbu.trace

From 153f490e6efbb5cd0ca1af5c89c01b22e410c8b9 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Tue, 24 Sep 2013 23:23:50 -0700
Subject: [PATCH 11/23] initial draft (spectrum similar to fbu preprint)

---
 python/generate_bump.py | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)
 create mode 100755 python/generate_bump.py

diff --git a/python/generate_bump.py b/python/generate_bump.py
new file mode 100755
index 0000000..3d57611
--- /dev/null
+++ b/python/generate_bump.py
@@ -0,0 +1,23 @@
+#!/bin/env python
+
+import numpy as np
+import matplotlib.mlab as mlab
+import matplotlib.pyplot as plt
+from math import sqrt
+
+offset, scale = 100.0, 500.0
+num_bins = 50
+x = np.random.exponential(scale=1.0, size=1e6) * scale + offset
+xm = np.random.normal(loc=4.0, scale=0.25, size=1e5) * scale + offset
+xtot = np.concatenate((x, xm))
+
+def smear(x, a=0.5, b=0.1) :
+   sigma = a*sqrt(x) + b*x
+   return x+np.random.normal(loc=0.0, scale=sigma)
+
+xsmear = [smear(x) for x in xtot]
+n, bins, patches = plt.hist(xtot, num_bins, facecolor='green', alpha=0.25, log=True)
+n, bins, patches = plt.hist(xsmear, num_bins, facecolor='blue', alpha=0.25, log=True)
+
+plt.show()
+plt.savefig('fallingExp_withBump.png')

From ab41ab7c039f4efba742b3c80a6cd19c2ca11d9e Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Mon, 7 Oct 2013 00:59:25 -0700
Subject: [PATCH 12/23] significant re-write

Details:
- factorize the code in well-defined functions
- add docstrings explaining what I am doing
- rewrite most of this script. Now plotting and using the resulting
histograms as bins/counts that can be used for the unfolding.
- compute normalized response matrix
- add plotting functions
---
 python/generate_bump.py | 61 +++++++++++++++++++++++++++++++++--------
 1 file changed, 49 insertions(+), 12 deletions(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 3d57611..892044c 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -5,19 +5,56 @@
 import matplotlib.pyplot as plt
 from math import sqrt
 
-offset, scale = 100.0, 500.0
-num_bins = 50
-x = np.random.exponential(scale=1.0, size=1e6) * scale + offset
-xm = np.random.normal(loc=4.0, scale=0.25, size=1e5) * scale + offset
-xtot = np.concatenate((x, xm))
+def generateTruthSpectrum(offset=100.0, scale=500.0, nbkg=1e6, nsig=1e5) :
+    """Generate pseudo-events with an observable in [offset, scale].
+    The observable spectrum is inspired to the one expected for a
+    di-jet resonance, where the background is distributed as a falling
+    exponential, and the signal is a gaussian on top of it.
+    """
+    bkg = np.random.exponential(scale=1.0, size=nbkg) * scale + offset
+    sig = np.random.normal(loc=4.0, scale=0.25, size=nsig) * scale + offset
+    tot = np.concatenate((bkg, sig))
+    return tot
 
 def smear(x, a=0.5, b=0.1) :
-   sigma = a*sqrt(x) + b*x
-   return x+np.random.normal(loc=0.0, scale=sigma)
+    """Smear the pseudo-events generated with generateTruthSpectrum as
+    to mimic the jet resolution, i.e. with a gaussianely-distributed
+    \delta, see eq. 13, eq. 14, and par. 7.1 of the FBU paper.
+    """
+    sigma = a*sqrt(x) + b*x
+    return x+np.random.normal(loc=0.0, scale=sigma)
 
-xsmear = [smear(x) for x in xtot]
-n, bins, patches = plt.hist(xtot, num_bins, facecolor='green', alpha=0.25, log=True)
-n, bins, patches = plt.hist(xsmear, num_bins, facecolor='blue', alpha=0.25, log=True)
+sigPlusBkg = generateTruthSpectrum()
+spbSmeared = [smear(x) for x in sigPlusBkg]
+numBins = 50
 
-plt.show()
-plt.savefig('fallingExp_withBump.png')
+nBinsTruth, truthMin, truthMax = 50, 500.0, 4500.0
+nBinsReco,  recoMin,  recoMax  = 50, 500.0, 4500.0
+histTruth, binsTruth        = np.histogram(sigPlusBkg, bins=nBinsTruth, range=(truthMin, truthMax))
+histReco,  binsReco         = np.histogram(spbSmeared, bins=nBinsReco, range=(recoMin, recoMax))
+respHist,  xedges, yedges   = np.histogram2d(sigPlusBkg, spbSmeared, bins=(nBinsTruth, nBinsReco),
+                                             range=((truthMin, truthMax), (recoMin,  recoMax)))
+def normalized(mat) :
+    mat = mat.astype(float)
+    return mat / mat.sum()
+respMat = normalized(respHist)
+print respHist
+print respMat
+
+def plotTruthAndReco(histT, binsT, histR, binsR, outfname) :
+    def getWidth(bins) : return bins[1] - bins[0]
+    plt.figure()
+    plt.bar(binsT[:-1], histT, color='b', width=getWidth(binsT), alpha=0.15)
+    plt.bar(binsR[:-1], histR, color='r', width=getWidth(binsR), alpha=0.15)
+    plt.savefig(outfname)
+
+def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
+    from matplotlib.colors import LogNorm
+    extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
+    plt.figure()
+    plt.imshow(respMatrix, extent=extent, interpolation='nearest', norm=LogNorm(), origin='lower')
+    plt.colorbar()
+    plt.savefig(outfname)
+
+plotTruthAndReco(histTruth, binsTruth, histReco, binsReco, 'fallingExp_withBump.png')
+plotRespMatrix(respHist, xedges, yedges, 'responseMatrix.png')

From 168ea4e0d4ac42cc61dd05f5582238852deea1b8 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Mon, 7 Oct 2013 01:21:10 -0700
Subject: [PATCH 13/23] utils to write/read np.array from/to json

---
 python/utils.py | 10 ++++++++++
 1 file changed, 10 insertions(+)
 create mode 100644 python/utils.py

diff --git a/python/utils.py b/python/utils.py
new file mode 100644
index 0000000..e409926
--- /dev/null
+++ b/python/utils.py
@@ -0,0 +1,10 @@
+import json
+import numpy as np
+
+def array2json(a, outfname) :
+    with open(outfname, 'w') as out:
+        json.dump(a.tolist(), out)
+
+def json2array(infname) :
+    with open(infname) as inp:
+        return numpy.array(json.load(inp))

From 3873eed813230a0d205e4dede057a7cc22bc4dc9 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Mon, 7 Oct 2013 01:42:11 -0700
Subject: [PATCH 14/23] separate the generation from the unfolding

Details:
two separate logic steps: generate the pseudo-data, and then unfold.
The pseudo-data can come either from the json file, or be generated.
---
 python/generate_bump.py | 80 +++++++++++++++++++++++++----------------
 1 file changed, 50 insertions(+), 30 deletions(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 892044c..2b4c6b9 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -4,6 +4,11 @@
 import matplotlib.mlab as mlab
 import matplotlib.pyplot as plt
 from math import sqrt
+import os
+from utils import array2json, json2array
+
+
+regenerate = True
 
 def generateTruthSpectrum(offset=100.0, scale=500.0, nbkg=1e6, nsig=1e5) :
     """Generate pseudo-events with an observable in [offset, scale].
@@ -24,37 +29,52 @@ def smear(x, a=0.5, b=0.1) :
     sigma = a*sqrt(x) + b*x
     return x+np.random.normal(loc=0.0, scale=sigma)
 
-sigPlusBkg = generateTruthSpectrum()
-spbSmeared = [smear(x) for x in sigPlusBkg]
-numBins = 50
+jsonDir = 'data/bump/'
+jsonTruthFname  = jsonDir+'truth.json'
+jsonSmearFname  = jsonDir+'reco.json'
+jsonResMatFname = jsonDir+'resMat.json'
 
-nBinsTruth, truthMin, truthMax = 50, 500.0, 4500.0
-nBinsReco,  recoMin,  recoMax  = 50, 500.0, 4500.0
-histTruth, binsTruth        = np.histogram(sigPlusBkg, bins=nBinsTruth, range=(truthMin, truthMax))
-histReco,  binsReco         = np.histogram(spbSmeared, bins=nBinsReco, range=(recoMin, recoMax))
-respHist,  xedges, yedges   = np.histogram2d(sigPlusBkg, spbSmeared, bins=(nBinsTruth, nBinsReco),
-                                             range=((truthMin, truthMax), (recoMin,  recoMax)))
 def normalized(mat) :
     mat = mat.astype(float)
     return mat / mat.sum()
-respMat = normalized(respHist)
-print respHist
-print respMat
-
-def plotTruthAndReco(histT, binsT, histR, binsR, outfname) :
-    def getWidth(bins) : return bins[1] - bins[0]
-    plt.figure()
-    plt.bar(binsT[:-1], histT, color='b', width=getWidth(binsT), alpha=0.15)
-    plt.bar(binsR[:-1], histR, color='r', width=getWidth(binsR), alpha=0.15)
-    plt.savefig(outfname)
-
-def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
-    from matplotlib.colors import LogNorm
-    extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
-    plt.figure()
-    plt.imshow(respMatrix, extent=extent, interpolation='nearest', norm=LogNorm(), origin='lower')
-    plt.colorbar()
-    plt.savefig(outfname)
-
-plotTruthAndReco(histTruth, binsTruth, histReco, binsReco, 'fallingExp_withBump.png')
-plotRespMatrix(respHist, xedges, yedges, 'responseMatrix.png')
+
+def generateAndPlot() :
+    sigPlusBkg = generateTruthSpectrum()
+    spbSmeared = np.array([smear(x) for x in sigPlusBkg])
+    numBins = 50
+    nBinsTruth, truthMin, truthMax = 50, 500.0, 4500.0
+    nBinsReco,  recoMin,  recoMax  = 50, 500.0, 4500.0
+    histTruth, binsTruth        = np.histogram(sigPlusBkg, bins=nBinsTruth, range=(truthMin, truthMax))
+    histReco,  binsReco         = np.histogram(spbSmeared, bins=nBinsReco, range=(recoMin, recoMax))
+    respHist,  xedges, yedges   = np.histogram2d(sigPlusBkg, spbSmeared, bins=(nBinsTruth, nBinsReco),
+                                                 range=((truthMin, truthMax), (recoMin,  recoMax)))
+    respMat = normalized(respHist)
+    def plotTruthAndReco(histT, binsT, histR, binsR, outfname) :
+        def getWidth(bins) : return bins[1] - bins[0]
+        plt.figure()
+        plt.bar(binsT[:-1], histT, color='b', width=getWidth(binsT), alpha=0.15)
+        plt.bar(binsR[:-1], histR, color='r', width=getWidth(binsR), alpha=0.15)
+        plt.savefig(outfname)
+    def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
+        from matplotlib.colors import LogNorm
+        extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
+        plt.figure()
+        plt.imshow(respMatrix, extent=extent, interpolation='nearest', norm=LogNorm(), origin='lower')
+        plt.colorbar()
+        plt.savefig(outfname)
+    plotTruthAndReco(histTruth, binsTruth, histReco, binsReco, 'fallingExp_withBump.png')
+    plotRespMatrix(respHist, xedges, yedges, 'responseMatrix.png')
+    return histTruth, histReco, respMat
+
+histTruth, histReco, respMat = None, None, None
+if regenerate :
+    histTruth, histReco, respMat = generateAndPlot()
+    if not os.path.isdir(jsonDir) : os.path.mkdir(jsonDir)
+    for a,f in [(histTruth, jsonTruthFname),
+                (histReco,  jsonSmearFname),
+                (respMat,   jsonResMatFname)] :
+        array2json(a,f)
+else :
+    histTruth  = json2array(jsonTruthFname)
+    histReco   = json2array(jsonSmearFname),
+    respMat    = json2array(jsonResMatFname)

From 5a5dfab8fb74d2dd77606fa36929d5b3c84c42c2 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Mon, 7 Oct 2013 01:46:46 -0700
Subject: [PATCH 15/23] almost there, but pyfbu does not work

Details:
complaining about template. I thought it wasn't needed anymore.
Need to ask clement and Francesco about it.
---
 python/generate_bump.py | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 2b4c6b9..18d0b47 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -6,7 +6,7 @@
 from math import sqrt
 import os
 from utils import array2json, json2array
-
+from pyFBU import pyFBU
 
 regenerate = True
 
@@ -78,3 +78,20 @@ def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
     histTruth  = json2array(jsonTruthFname)
     histReco   = json2array(jsonSmearFname),
     respMat    = json2array(jsonResMatFname)
+
+pyfbu = pyFBU()
+pyfbu.nMCMC        = 100000
+pyfbu.nBurn        = 1000
+pyfbu.nThin        = 10
+pyfbu.lower        = 70000
+pyfbu.upper        = 140000
+pyfbu.jsonData     = jsonSmearFname
+pyfbu.jsonMig      = jsonResMatFname
+pyfbu.jsonBkg      = '' # this shouldnt be a requirement; pass [0.] ?
+pyfbu.templateFile = '' # templateFile  ?
+pyfbu.modelName    = '' # modelName if modelName else pyfbu.modelName ? I thought this wasn't necessary anymore
+pyfbu.verbose      = True # verbose ? will be on cmd-line
+
+pyfbu.run()
+
+trace = pyfbu.trace

From 580ed110d78e6ab7367559749462bcdf8fcdd662 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Wed, 9 Oct 2013 21:19:25 -0700
Subject: [PATCH 16/23] flip up-down the response matrix

Details:
the vertical axis of the histo provided by histogram2d is reversed
(see
http://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram2d.html)
---
 python/generate_bump.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 18d0b47..0f81787 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -49,6 +49,7 @@ def generateAndPlot() :
     respHist,  xedges, yedges   = np.histogram2d(sigPlusBkg, spbSmeared, bins=(nBinsTruth, nBinsReco),
                                                  range=((truthMin, truthMax), (recoMin,  recoMax)))
     respMat = normalized(respHist)
+    respMat = np.flipud(respMat) # histogram2d returns the y axis reversed (see docs)
     def plotTruthAndReco(histT, binsT, histR, binsR, outfname) :
         def getWidth(bins) : return bins[1] - bins[0]
         plt.figure()

From 001ecfff0ec64f7a0504b1e127e691497ba68bdd Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Wed, 9 Oct 2013 21:21:41 -0700
Subject: [PATCH 17/23] delete assignement of obsolete parameters

---
 python/generate_bump.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 0f81787..5aac03d 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -88,9 +88,6 @@ def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
 pyfbu.upper        = 140000
 pyfbu.jsonData     = jsonSmearFname
 pyfbu.jsonMig      = jsonResMatFname
-pyfbu.jsonBkg      = '' # this shouldnt be a requirement; pass [0.] ?
-pyfbu.templateFile = '' # templateFile  ?
-pyfbu.modelName    = '' # modelName if modelName else pyfbu.modelName ? I thought this wasn't necessary anymore
 pyfbu.verbose      = True # verbose ? will be on cmd-line
 
 pyfbu.run()

From 9710442f94ba3937b50fa2e5ffe9f764bbd5f4da Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Wed, 9 Oct 2013 21:22:27 -0700
Subject: [PATCH 18/23] reduce the number of samplings to test things out
 faster

---
 python/generate_bump.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 5aac03d..102779f 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -81,8 +81,8 @@ def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
     respMat    = json2array(jsonResMatFname)
 
 pyfbu = pyFBU()
-pyfbu.nMCMC        = 100000
-pyfbu.nBurn        = 1000
+pyfbu.nMCMC        = 1000
+pyfbu.nBurn        = 100
 pyfbu.nThin        = 10
 pyfbu.lower        = 70000
 pyfbu.upper        = 140000

From 4106129ca82d0604a8626e3426525f03e94f81f6 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Wed, 9 Oct 2013 21:28:10 -0700
Subject: [PATCH 19/23] typo

---
 python/utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/utils.py b/python/utils.py
index e409926..61cbd43 100644
--- a/python/utils.py
+++ b/python/utils.py
@@ -7,4 +7,4 @@ def array2json(a, outfname) :
 
 def json2array(infname) :
     with open(infname) as inp:
-        return numpy.array(json.load(inp))
+        return np.array(json.load(inp))

From 69b0b6bcc94ceeed7e48b1e719ba851f8a528f7c Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Wed, 9 Oct 2013 21:34:14 -0700
Subject: [PATCH 20/23] regenerate when json is not there

---
 python/generate_bump.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 102779f..1a70472 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -8,8 +8,6 @@
 from utils import array2json, json2array
 from pyFBU import pyFBU
 
-regenerate = True
-
 def generateTruthSpectrum(offset=100.0, scale=500.0, nbkg=1e6, nsig=1e5) :
     """Generate pseudo-events with an observable in [offset, scale].
     The observable spectrum is inspired to the one expected for a
@@ -34,6 +32,8 @@ def smear(x, a=0.5, b=0.1) :
 jsonSmearFname  = jsonDir+'reco.json'
 jsonResMatFname = jsonDir+'resMat.json'
 
+regenerate = not os.path.exists(jsonResMatFname)
+
 def normalized(mat) :
     mat = mat.astype(float)
     return mat / mat.sum()

From ab547b222594eefa66ae781fb5de25224c1add06 Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Wed, 9 Oct 2013 21:34:41 -0700
Subject: [PATCH 21/23] try to adjust lower/upper

---
 python/generate_bump.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 1a70472..33e3a53 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -84,7 +84,7 @@ def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
 pyfbu.nMCMC        = 1000
 pyfbu.nBurn        = 100
 pyfbu.nThin        = 10
-pyfbu.lower        = 70000
+pyfbu.lower        = 1
 pyfbu.upper        = 140000
 pyfbu.jsonData     = jsonSmearFname
 pyfbu.jsonMig      = jsonResMatFname

From dbe9c3157aa2db4a2795cbe0be9eabcd4f30a11d Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Sun, 13 Oct 2013 04:38:10 -0700
Subject: [PATCH 22/23] still work in progress, getting there

---
 python/generate_bump.py | 57 ++++++++++++++++++++++++++++++++++++-----
 1 file changed, 51 insertions(+), 6 deletions(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index 33e3a53..e4bea56 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -8,7 +8,11 @@
 from utils import array2json, json2array
 from pyFBU import pyFBU
 
+nBinsTruth, truthMin, truthMax = 12, 500.0, 4500.0
+nBinsReco,  recoMin,  recoMax  = 12, 500.0, 4500.0
+
 def generateTruthSpectrum(offset=100.0, scale=500.0, nbkg=1e6, nsig=1e5) :
+# def generateTruthSpectrum(offset=100.0, scale=500.0, nbkg=10, nsig=2) :
     """Generate pseudo-events with an observable in [offset, scale].
     The observable spectrum is inspired to the one expected for a
     di-jet resonance, where the background is distributed as a falling
@@ -33,6 +37,7 @@ def smear(x, a=0.5, b=0.1) :
 jsonResMatFname = jsonDir+'resMat.json'
 
 regenerate = not os.path.exists(jsonResMatFname)
+regenerate = True
 
 def normalized(mat) :
     mat = mat.astype(float)
@@ -41,15 +46,36 @@ def normalized(mat) :
 def generateAndPlot() :
     sigPlusBkg = generateTruthSpectrum()
     spbSmeared = np.array([smear(x) for x in sigPlusBkg])
-    numBins = 50
-    nBinsTruth, truthMin, truthMax = 50, 500.0, 4500.0
-    nBinsReco,  recoMin,  recoMax  = 50, 500.0, 4500.0
     histTruth, binsTruth        = np.histogram(sigPlusBkg, bins=nBinsTruth, range=(truthMin, truthMax))
     histReco,  binsReco         = np.histogram(spbSmeared, bins=nBinsReco, range=(recoMin, recoMax))
     respHist,  xedges, yedges   = np.histogram2d(sigPlusBkg, spbSmeared, bins=(nBinsTruth, nBinsReco),
                                                  range=((truthMin, truthMax), (recoMin,  recoMax)))
-    respMat = normalized(respHist)
-    respMat = np.flipud(respMat) # histogram2d returns the y axis reversed (see docs)
+    print 'in truth histo ',histTruth.shape,': ',histTruth
+    print 'in reco histo  ',histReco.shape ,': ',histReco
+    print 'in respHist    ',respHist.shape ,': ',respHist
+    def prMatrix(recoVsTruthCounts) :
+        """
+        This matrix is p(r|t). Given recoVsTruthCounts, where truth is
+        on axis=0 and reco is on axis=1, we can just obtain the p(r|t)
+        matrix by normalizing each bin [t][r] by the total N of evts
+        in that 'row' [t]:
+
+           [t][r] -> [t][r]/sum([[t][r] for r in nreco])
+
+        With np this can be done with the broadcasting mechanims. Try:
+        a = np.array([[1., 2.],
+                      [3., 4.]])
+        a / np.sum(a, axis=1)[:,np.newaxis]
+        --> [[ 0.33333333,  0.66666667],
+             [ 0.42857143,  0.57142857]]
+        And see for example
+        http://stackoverflow.com/questions/7140738/numpy-divide-along-axis)
+        """
+        return recoVsTruthCounts / np.sum(recoVsTruthCounts, axis=1)[:,np.newaxis]
+#return recoVsTruthCounts / truthCounts.astype('float')[:,np.newaxis]
+    respMat = prMatrix(respHist)
+    print 'respMat :'
+    print respMat
     def plotTruthAndReco(histT, binsT, histR, binsR, outfname) :
         def getWidth(bins) : return bins[1] - bins[0]
         plt.figure()
@@ -80,8 +106,11 @@ def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
     histReco   = json2array(jsonSmearFname),
     respMat    = json2array(jsonResMatFname)
 
+print 'histTruth : ',histTruth
+print 'histReco : ',histReco
+print 'respMat : ',respMat
 pyfbu = pyFBU()
-pyfbu.nMCMC        = 1000
+pyfbu.nMCMC        = 100000
 pyfbu.nBurn        = 100
 pyfbu.nThin        = 10
 pyfbu.lower        = 1
@@ -93,3 +122,19 @@ def plotRespMatrix(respMatrix, xedges, yedges, outfname) :
 pyfbu.run()
 
 trace = pyfbu.trace
+print trace
+estMean   = np.mean(trace, axis=0)
+estMedian = np.median(trace, axis=0)
+estStd    = np.std(trace, axis=0)
+print 'estMean: ',estMean
+print 'estMedian: ',estMedian
+xdata = np.linspace(truthMin, truthMax, nBinsTruth)
+binW = (truthMax-truthMin)/nBinsTruth
+plt.figure()
+print 'xdata ',len(xdata)
+print 'histTruth ',len(histTruth)
+plt.bar(xdata, histTruth, color='b', width=binW, alpha=0.15)
+#plt.bar(xdata, histReco, color='r',  width=binW, alpha=0.15)
+
+plt.errorbar(xdata+0.5*binW, estMedian, yerr=estStd, fmt='k.')
+plt.savefig('foo.png')

From 0587d84b7072cc39bedb7f5db9cb01e5b25ee67d Mon Sep 17 00:00:00 2001
From: Davide Gerbaudo <davide.gerbaudo@gmail.com>
Date: Mon, 18 Nov 2013 10:41:27 -0800
Subject: [PATCH 23/23] protect against divide-by-zero

Details:
when computing p(r|t), there can be bins without any reconstructed
events. Protect against divide by zero.
---
 python/generate_bump.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/python/generate_bump.py b/python/generate_bump.py
index e4bea56..48e3d4f 100755
--- a/python/generate_bump.py
+++ b/python/generate_bump.py
@@ -71,7 +71,9 @@ def prMatrix(recoVsTruthCounts) :
         And see for example
         http://stackoverflow.com/questions/7140738/numpy-divide-along-axis)
         """
-        return recoVsTruthCounts / np.sum(recoVsTruthCounts, axis=1)[:,np.newaxis]
+        den = np.sum(recoVsTruthCounts, axis=1)[:,np.newaxis]
+        den[den==0.0]=1.0 # avoid divide by zero; if all values in this row are zero, 0./1. is still 0.
+        return recoVsTruthCounts / den
 #return recoVsTruthCounts / truthCounts.astype('float')[:,np.newaxis]
     respMat = prMatrix(respHist)
     print 'respMat :'