diff --git a/examples/discriminationMC.py b/examples/discriminationMC.py
index b2a57ef..b36cedd 100644
--- a/examples/discriminationMC.py
+++ b/examples/discriminationMC.py
@@ -1,185 +1,194 @@
 '''
 discriminationMC.py
 Demonstration of the discriminating power of the debinning calculation.
 
 Copyright (C) 2015 Abram Krislock
 
 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program.  If not, see http://www.gnu.org/licenses/gpl-3.0.txt.
 '''
 
 ### Imports
 from __future__ import division
 from matplotlib import pyplot as plt
 import numpy as np
 from collections import OrderedDict as oDict
 import bisect
 import os
 
 debinningDir = os.getcwd()
 while debinningDir.rfind('/') > debinningDir.rfind('debinning'):
   debinningDir = debinningDir[:debinningDir.rfind('/')]
 if debinningDir not in os.sys.path:
   os.sys.path.insert(0, debinningDir)
 
 from debinningCore import settings
 settings.useGPU = True
 settings.invertCovariance = True
 import debinningCore.debinning as debinning
 import debinningMath.orderStatistics as oStat
 import debinningMath.dataGen as dataGen
 import debinningMath.sampleFunctions as funcs
 
 def iseDiscrimination():
   fig, ((ax1,ax2),(ax3,ax4)) = plt.subplots(nrows=2, ncols=2, figsize=(22,14),
                                             facecolor='#ffffff')
   bggen = dataGen.DataGenContainer(funcs.endpointBG,
                                    domain=np.linspace(0., 200., 256))
   mypoisson = dataGen.MyRandom('mypoisson')
   dataGen.randomGenerator.setSeed('someseed')
 
   nExpected = [100, 200, 500, 1000]
+  axes = {100:ax1, 200:ax2, 500:ax3, 1000:ax4}
   alpha = oDict((('zero',0.), ('one',0.01), ('three',0.03), ('seven',0.07)))
 # alpha = oDict((('zero',0.), ('one',0.01), ('two',0.02), ('three',0.03),
 #                ('four',0.04), ('five',0.05), ('six',0.06), ('seven',0.07),
 #                ('eight',0.08), ('nine',0.09)))
   debinChisq = {}
-  debinChisqMod = {}
   debinCovarChisq = {}
+  debinChisqMod = {}
   debinCovarChisqMod = {}
   neg2Loglike = {}
 
   for nEx in nExpected:
     for key in alpha:
       print 'For nExpected = ' + str(nEx) + ', alpha = ' + key + ' %'
       gen = dataGen.DataPlusBGContainer(funcs.endpoint, funcs.endpointBG,
                                         alpha[key], nPoints=nEx,
                                         domain=np.linspace(0., 200., 256))
       debinCalc = debinning.debinningCalculator(gen, bggen, nEx)
       dx = debinCalc.fit.x[1] - debinCalc.fit.x[0]
 
       debinChisq[key + str(nEx)] = []
-      debinChisqMod[key + str(nEx)] = []
       debinCovarChisq[key + str(nEx)] = []
+      debinChisqMod[key + str(nEx)] = []
       debinCovarChisqMod[key + str(nEx)] = []
       neg2Loglike[key + str(nEx)] = []
 
-      for n in mypoisson.poisson(nEx / (1. - alpha[key]), 10000):
+      for i in xrange(10000):
+        n = -1
+        while n < nEx:
+          n = mypoisson.poisson(nEx / (1. - alpha[key]))
         # NOTE: "True" n for data generation, nEx for debinning calculation
         debinCalc.gen.n = n
         debinCalc.resetData()
         debinCalc.gen.n = nEx
         debinCalc.calculateDebinning()
         # NOTE: The bg with nEx events is assumed to be a "good fit"
         debinCalc.calculateExpectations(True)  # <= send in True for goodFit
 
-        difference = (debinCalc.debinData - debinCalc.fit.f)
-        differMod = (n*debinCalc.debinData - nEx*debinCalc.fit.f)
+        difference = debinCalc.debinData - debinCalc.fit.f
+        differMod = n*debinCalc.debinData - nEx*debinCalc.fit.f
         debinChisq[key + str(nEx)].append(
                 np.sum(difference**2 / debinCalc.sigmaSq_x)
                 -2.* (n * np.log(nEx) - nEx - oStat.logfact(n)) )
-        debinChisqMod[key + str(nEx)].append(
-                np.sum(differMod**2 / debinCalc.sigmaSq_x / nEx))
         debinCovarChisq[key + str(nEx)].append(
                 np.dot(difference, np.dot(debinCalc.invV_xy, difference))
                 -2.* (n * np.log(nEx) - nEx - oStat.logfact(n)) )
+        debinChisqMod[key + str(nEx)].append(
+                np.sum(differMod**2 / debinCalc.sigmaSq_x) / nEx )
         debinCovarChisqMod[key + str(nEx)].append(
-                np.dot(differMod, np.dot(debinCalc.invV_xy, differMod)) / nEx)
+                np.dot(differMod, np.dot(debinCalc.invV_xy, differMod)) / nEx )
         neg2Loglike[key + str(nEx)].append(
                 -2. * (np.log(debinCalc.fit.f.take(debinCalc.dataToX)).sum()
                 + n * np.log(nEx) - nEx - oStat.logfact(n)) )
 
       # Calculate empirical distribution functions for power calculation
       neg2Loglike[key + str(nEx)] = np.sort(neg2Loglike[key + str(nEx)])
       debinChisq[key + str(nEx)] = np.sort(debinChisq[key + str(nEx)])
-      debinChisqMod[key + str(nEx)] = np.sort(debinChisqMod[key + str(nEx)])
       debinCovarChisq[key + str(nEx)] = np.sort(debinCovarChisq[key + str(nEx)])
+      debinChisqMod[key + str(nEx)] = np.sort(debinChisqMod[key + str(nEx)])
       debinCovarChisqMod[key + str(nEx)] = np.sort(debinCovarChisqMod[key + str(nEx)])
 
       # Power and p-value calculations
       if key == 'zero':
         # Once the data is sorted, the 95th percentile corresponds to the 95th 
         #   percentile of the indices. 10000 * 0.95 = 9500 (minus one).
         e95Loglike = neg2Loglike[key + str(nEx)][9499]
         e95Chisq = debinChisq[key + str(nEx)][9499]
-        e95ChisqMod = debinChisqMod[key + str(nEx)][9499]
         e95CovarChisq = debinCovarChisq[key + str(nEx)][9499]
+        e95ChisqMod = debinChisqMod[key + str(nEx)][9499]
         e95CovarChisqMod = debinCovarChisqMod[key + str(nEx)][9499]
       else:
         try:
           # bisect_right(data, m) locates the right-most entry in data which is <= m
           beta = bisect.bisect_right(neg2Loglike[key + str(nEx)], e95Loglike) / 10000
           print '  -2 log like power at 95% exclusion = ' + str(1. - beta)
         except IndexError:
           print '  WTF for -2 log like: ' + str(e95Loglike)
           print
           print neg2Loglike[key + str(nEx)]
           print
           print
         try:
           beta = bisect.bisect_right(debinChisq[key + str(nEx)], e95Chisq) / 10000
           print '  debinning chisq power at 95% exclusion = ' + str(1. - beta)
         except IndexError:
           print '  WTF for debinning: ' + str(e95Chisq)
           print
           print debinChisq[key + str(nEx)]
           print
           print
         try:
-          beta = bisect.bisect_right(debinChisqMod[key + str(nEx)], e95ChisqMod) / 10000
-          print '  debinning chisqMod power at 95% exclusion = ' + str(1. - beta)
+          beta = bisect.bisect_right(debinCovarChisq[key + str(nEx)], e95CovarChisq) / 10000
+          print '  debinning CovarChisq power at 95% exclusion = ' + str(1. - beta)
         except IndexError:
-          print '  WTF for debinning: ' + str(e95ChisqMod)
+          print '  WTF for debinning: ' + str(e95CovarChisq)
           print
-          print debinChisqMod[key + str(nEx)]
+          print debinCovarChisq[key + str(nEx)]
           print
           print
         try:
-          beta = bisect.bisect_right(debinCovarChisq[key + str(nEx)], e95CovarChisq) / 10000
-          print '  debinning CovarChisq power at 95% exclusion = ' + str(1. - beta)
+          beta = bisect.bisect_right(debinChisqMod[key + str(nEx)], e95ChisqMod) / 10000
+          print '  debinning chisqMod power at 95% exclusion = ' + str(1. - beta)
         except IndexError:
-          print '  WTF for debinning: ' + str(e95CovarChisq)
+          print '  WTF for debinning: ' + str(e95ChisqMod)
           print
-          print debinCovarChisq[key + str(nEx)]
+          print debinChisqMod[key + str(nEx)]
           print
           print
         try:
           beta = bisect.bisect_right(debinCovarChisqMod[key + str(nEx)], e95CovarChisqMod) / 10000
           print '  debinning CovarChisqMod power at 95% exclusion = ' + str(1. - beta)
         except IndexError:
           print '  WTF for debinning: ' + str(e95CovarChisqMod)
           print
           print debinCovarChisqMod[key + str(nEx)]
           print
           print
       print
 
+      axes[nEx].hist(debinCovarChisq[key + str(nEx)], bins=51, alpha=0.4,
+               label=r'$\alpha = %f$' % alpha[key])
+      axes[nEx].set_xlabel('Covariant Chisq, nEx=' + str(nEx), labelpad=5, fontsize=22)
+'''
       if nEx == 1000:
         ax1.hist(debinChisq[key + str(nEx)], bins=51, alpha=0.4,
                  label=r'$\alpha = %f$' % alpha[key])
         ax1.set_xlabel('Chisq, nEx=1000', labelpad=5, fontsize=22)
-        ax2.hist(debinChisqMod[key + str(nEx)], bins=51, alpha=0.4,
+        ax2.hist(debinCovarChisq[key + str(nEx)], bins=51, alpha=0.4,
                  label=r'$\alpha = %f$' % alpha[key])
-        ax2.set_xlabel('ChisqMod, nEx=1000', labelpad=5, fontsize=22)
-        ax3.hist(debinCovarChisq[key + str(nEx)], bins=51, alpha=0.4,
+        ax2.set_xlabel('CovarChisq, nEx=1000', labelpad=5, fontsize=22)
+        ax3.hist(debinChisqMod[key + str(nEx)], bins=51, alpha=0.4,
                  label=r'$\alpha = %f$' % alpha[key])
-        ax3.set_xlabel('CovarChisq, nEx=1000', labelpad=5, fontsize=22)
+        ax3.set_xlabel('ChisqMod, nEx=1000', labelpad=5, fontsize=22)
         ax4.hist(debinCovarChisqMod[key + str(nEx)], bins=51, alpha=0.4,
                  label=r'$\alpha = %f$' % alpha[key])
         ax4.set_xlabel('CovarChisqMod, nEx=1000', labelpad=5, fontsize=22)
+'''
       
 if __name__ == '__main__':
   print 'Now running the discrimination demonstration for the endpointBG function'
   plt.ion()
   iseDiscrimination()
   print 'Done.\n\n'