diff options
Diffstat (limited to 'utils/chi2')
| -rwxr-xr-x | utils/chi2 | 16 |
1 files changed, 9 insertions, 7 deletions
@@ -189,7 +189,7 @@ FIT_PARS = [ 'Electron + Muon energy resolution', 'External Muon scale'] -def make_nll(data, muons, mc, bins): +def make_nll(data, muons, mc, bins, print_nll=False): data_hists = get_data_hists(data,bins) ke_dict = {} @@ -257,8 +257,9 @@ def make_nll(data, muons, mc, bins): nll -= norm.logpdf(x[4],ENERGY_RESOLUTION_MEAN['u'],ENERGY_RESOLUTION_UNCERTAINTY['u']) nll -= norm.logpdf(x[6],ENERGY_RESOLUTION_MEAN['eu'],ENERGY_RESOLUTION_UNCERTAINTY['eu']) - # Print the result - print("nll = %.2f" % nll) + if print_nll: + # Print the result + print("nll = %.2f" % nll) return nll return nll @@ -359,7 +360,7 @@ def get_prob(data,muon,mc,samples,bins,size): print(id, prob[id]) return prob -def do_fit(data,muon,data_mc,bins,steps): +def do_fit(data,muon,data_mc,bins,steps,print_nll=False): """ Run the fit and return the minimum along with samples from running an MCMC starting near the minimum. @@ -376,7 +377,7 @@ def do_fit(data,muon,data_mc,bins,steps): Returns a tuple (xopt, samples) where samples is an array of shape (steps, number of parameters). """ - nll = make_nll(data,muon,data_mc,bins) + nll = make_nll(data,muon,data_mc,bins,print_nll) x0 = np.array([1.0,1.0,EPSILON,1.0,EPSILON,1.0,EPSILON,EPSILON]) opt = nlopt.opt(nlopt.LN_SBPLX, len(x0)) @@ -442,6 +443,7 @@ if __name__ == '__main__': parser.add_argument("--multinomial-prob-size", type=int, default=10000, help="number of p values to compute") parser.add_argument("--coverage", type=int, default=0, help="plot p value coverage") parser.add_argument("--weights", nargs='+', required=True, help="GENIE reweight HDF5 files") + parser.add_argument("--print-nll", action='store_true', default=False, help="print nll values") args = parser.parse_args() setup_matplotlib(args.save) @@ -558,7 +560,7 @@ if __name__ == '__main__': data.ke += np.random.randn(len(data.ke))*data.ke*energy_resolution data_atm.ke += np.random.randn(len(data_atm.ke))*data_atm.ke*energy_resolution - xopt, samples = do_fit(data,muon,data_mc,bins,args.steps) + xopt, samples = do_fit(data,muon,data_mc,bins,args.steps,args.print_nll) for i in range(len(FIT_PARS)): mean = np.mean(samples[:,i]) @@ -643,7 +645,7 @@ if __name__ == '__main__': sys.exit(0) - xopt, samples = do_fit(data,muon,data_mc,bins,args.steps) + xopt, samples = do_fit(data,muon,data_mc,bins,args.steps,args.print_nll) prob = get_prob(data,muon,data_mc_with_weights,samples,bins,size=args.multinomial_prob_size) prob_atm = get_prob(data_atm,muon_atm,data_atm_mc_with_weights,samples,bins,size=args.multinomial_prob_size) |