update contamination analysis stuff

- fix Constraint.renormalize_no_fix() which could enter an infinite loop if the
  fixed parameter was greater than 1 - EPSILON
- don't divide by psi twice in get_events()
- only use prompt events and cut on nhit_cal < 100

4 files changed, 31 insertions, 9 deletions

diff --git a/utils/dc b/utils/dc
index 834c6e5..8227e99 100755
--- a/utils/dc
+++ b/utils/dc
@@ -280,6 +280,8 @@ if __name__ == '__main__':
     import matplotlib.pyplot as plt
 
     ev = get_events(args.filenames,merge_fits=True)
+    ev = ev[ev.prompt]
+    ev = ev[ev.nhit_cal > 100]
 
     # figure out bins for high level variables
     ev = radius_cut(ev)
@@ -451,7 +453,7 @@ if __name__ == '__main__':
     #samples = metropolis_hastings(nll,xopt,stepsizes,100000)
     #print("nll(xopt) = %.2g" % nll(xopt))
 
-    pos = np.empty((10, len(x0)),dtype=np.double)
+    pos = np.empty((100, len(x0)),dtype=np.double)
     for i in range(pos.shape[0]):
         pos[i] = xopt + np.random.randn(len(x0))*stepsizes
         pos[i,:6] = np.clip(pos[i,:6],EPSILON,1e9)
@@ -463,7 +465,7 @@ if __name__ == '__main__':
 
     nwalkers, ndim = pos.shape
 
-    proposal = get_proposal_func(stepsizes*0.1,low,high)
+    proposal = get_proposal_func(stepsizes,low,high)
     sampler = emcee.EnsembleSampler(nwalkers, ndim, lambda x, grad, fill_value: -nll(x,grad,fill_value), moves=emcee.moves.MHMove(proposal),args=[None,np.inf])
     with np.errstate(invalid='ignore'):
         sampler.run_mcmc(pos, args.steps)
diff --git a/utils/dc-closure-test b/utils/dc-closure-test
index 54d2629..45c0f6e 100755
--- a/utils/dc-closure-test
+++ b/utils/dc-closure-test
@@ -383,7 +383,7 @@ def fit(data, sacrifice, steps):
     #samples = metropolis_hastings(nll,xopt,stepsizes,100000)
     #print("nll(xopt) = %.2g" % nll(xopt))
 
-    pos = np.empty((10, len(x0)),dtype=np.double)
+    pos = np.empty((100, len(x0)),dtype=np.double)
     for i in range(pos.shape[0]):
         pos[i] = xopt + np.random.randn(len(x0))*stepsizes
         pos[i,:6] = np.clip(pos[i,:6],EPSILON,1e9)
@@ -395,7 +395,7 @@ def fit(data, sacrifice, steps):
 
     nwalkers, ndim = pos.shape
 
-    proposal = get_proposal_func(stepsizes*0.1,low,high)
+    proposal = get_proposal_func(stepsizes,low,high)
     sampler = emcee.EnsembleSampler(nwalkers, ndim, lambda x, grad, fill_value: -nll(x,grad,fill_value), moves=emcee.moves.MHMove(proposal),args=[None,np.inf])
     with np.errstate(invalid='ignore'):
         sampler.run_mcmc(pos, steps)
@@ -425,7 +425,8 @@ if __name__ == '__main__':
     import matplotlib.pyplot as plt
 
     ev = get_events(args.filenames,merge_fits=True)
-    ev_mc = get_events(args.mc, merge_fits=True)
+    ev = ev[ev.prompt]
+    ev = ev[ev.nhit_cal > 100]
 
     # figure out bins for high level variables
     ev = radius_cut(ev)
@@ -441,7 +442,9 @@ if __name__ == '__main__':
     ev['muon'] = ((ev.dc & DC_MUON) != 0) & ~(ev.noise | ev.neck | ev.flasher | ev.breakdown)
     ev['signal'] = ~(ev.noise | ev.neck | ev.flasher | ev.breakdown | ev.muon)
 
+    ev_mc = get_events(args.mc, merge_fits=True)
     ev_mc = ev_mc[ev_mc.prompt]
+    ev_mc = ev_mc[ev_mc.nhit_cal > 100]
 
     # figure out bins for high level variables
     ev_mc = radius_cut(ev_mc)
diff --git a/utils/sddm/dc.py b/utils/sddm/dc.py
index a1d87b5..d681831 100755
--- a/utils/sddm/dc.py
+++ b/utils/sddm/dc.py
@@ -343,12 +343,29 @@ class Constraint(object):
     def renormalize(self, x, fix):
         if fix not in self.range:
             return x
+
+        if 1 - x[fix] < EPSILON:
+            return x
+
+        if self(x) < 0:
+            return x
+
+        # Scale all the other entries except for the fixed one by a constant to
+        # ensure that sum(x) = 1 - EPSILON
+        c = ((1-x[fix])-EPSILON)/sum(x[i] for i in self.range if i != fix)
+
         x = x.copy()
+
+        for i in self.range:
+            if i != fix:
+                x[i] *= c
+
         while self(x) >= 0:
             for i in self.range:
                 if i == fix:
                     continue
-                x[i] -= self(x)
+                x[i] /= 2
+
         return x
 
     def renormalize_no_fix(self, x):
diff --git a/utils/sddm/plot_energy.py b/utils/sddm/plot_energy.py
index c9f5142..88950f7 100755
--- a/utils/sddm/plot_energy.py
+++ b/utils/sddm/plot_energy.py
@@ -471,13 +471,13 @@ def get_events(filenames, merge_fits=False):
     ev = ev.groupby('run',group_keys=False).apply(retrigger_cut)
 
     if merge_fits:
-        ev.set_index(['run','gtid'])
-
         ev = pd.merge(fits,ev,how='inner',on=['run','gtid'])
+        # Set the index to (run, gtid) so we can set columns from the single particle results
+        ev = ev.set_index(['run','gtid'])
         ev_single_particle = ev[(ev.id2 == 0) & (ev.id3 == 0)]
         ev_single_particle = ev_single_particle.sort_values('fmin').groupby(['run','gtid']).nth(0)
         ev = ev.sort_values('fmin').groupby(['run','gtid']).nth(0)
-        ev['psi'] /= ev.nhit_cal
+        #ev['psi'] /= ev.nhit_cal
 
         ev['cos_theta'] = np.cos(ev_single_particle['theta1'])
         ev['r'] = np.sqrt(ev.x**2 + ev.y**2 + ev.z**2)