diff options
Diffstat (limited to 'utils/plot-michels')
| -rwxr-xr-x | utils/plot-michels | 22 |
1 files changed, 9 insertions, 13 deletions
diff --git a/utils/plot-michels b/utils/plot-michels index b313c3e..9637e13 100755 --- a/utils/plot-michels +++ b/utils/plot-michels @@ -30,17 +30,7 @@ from sddm.stats import * import emcee from sddm.dc import estimate_errors, EPSILON import nlopt -import contextlib - -# from https://stackoverflow.com/questions/2891790/how-to-pretty-print-a-numpy-array-without-scientific-notation-and-with-given-pre -@contextlib.contextmanager -def printoptions(*args, **kwargs): - original = np.get_printoptions() - np.set_printoptions(*args, **kwargs) - try: - yield - finally: - np.set_printoptions(**original) +from sddm import printoptions particle_id = {20: 'e', 22: 'u'} @@ -100,6 +90,7 @@ if __name__ == '__main__': from sddm.plot_energy import * from sddm.plot import * from sddm import setup_matplotlib + from sddm.utils import correct_energy_bias parser = argparse.ArgumentParser("plot fit results") parser.add_argument("filenames", nargs='+', help="input files") @@ -108,6 +99,7 @@ if __name__ == '__main__': parser.add_argument("--nhit-thresh", type=int, default=None, help="nhit threshold to apply to events before processing (should only be used for testing to speed things up)") parser.add_argument("--print-nll", action='store_true', default=False, help="print nll values") parser.add_argument("--steps", type=int, default=1000, help="number of steps in the MCMC chain") + parser.add_argument("--walkers", type=int, default=100, help="number of walkers") args = parser.parse_args() setup_matplotlib(args.save) @@ -120,7 +112,9 @@ if __name__ == '__main__': for run, df in rhdr.groupby('run'): evs.append(get_events(df.filename.values, merge_fits=True, nhit_thresh=args.nhit_thresh)) ev = pd.concat(evs) + ev = correct_energy_bias(ev) ev_mc = get_events(args.mc, merge_fits=True) + ev_mc = correct_energy_bias(ev_mc) # Drop events without fits ev = ev[~np.isnan(ev.fmin)] @@ -195,17 +189,19 @@ if __name__ == '__main__': stepsizes = estimate_errors(nll,xopt,low,high) + stepsizes[1] = 0.1 + with printoptions(precision=3, suppress=True): print("Errors: ", stepsizes) - pos = np.empty((20, len(x0)),dtype=np.double) + pos = np.empty((args.walkers, len(x0)),dtype=np.double) for i in range(pos.shape[0]): pos[i] = xopt + np.random.randn(len(x0))*stepsizes pos[i,:] = np.clip(pos[i,:],low,high) nwalkers, ndim = pos.shape - sampler = emcee.EnsembleSampler(nwalkers, ndim, lambda x: -nll(x)) + sampler = emcee.EnsembleSampler(nwalkers, ndim, lambda x: -nll(x), moves=emcee.moves.KDEMove()) with np.errstate(invalid='ignore'): sampler.run_mcmc(pos, args.steps) |