From 604828af0f1f9411477976cf16648270ef1a0128 Mon Sep 17 00:00:00 2001 From: davidwalter2 Date: Thu, 19 Dec 2024 13:23:17 -0500 Subject: [PATCH] Latest updates on narf --- .github/workflows/main.yml | 2 +- .github/workflows/unfolding.yml | 38 ++++----- narf | 2 +- scripts/plotting/postfitPlots.py | 130 +++++++++++++++++++++++-------- 4 files changed, 119 insertions(+), 53 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 63e1bea79..b7befcf58 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -528,7 +528,7 @@ jobs: scripts/ci/run_with_singularity.sh scripts/ci/setup_and_run_python.sh narf/scripts/fitting/combinetf2.py $COMBINED_DIR/Combination_WMass_lowPUZMass_lowPU_ptVGen/Combination.hdf5 -o $COMBINED_DIR/Combination_WMass_lowPUZMass_lowPU_ptVGen/fitresults.hdf5 --externalPostfit $COMBINED_DIR/Combination_WMass_lowPUZMass_lowPU_ptll/fitresults.hdf5 - --saveHists --computeHistErrors + --saveHists --computeHistErrors --computeHistCov --computeHistImpacts --noChi2 - name: lowpu combine unfolded xsec plot run: >- diff --git a/.github/workflows/unfolding.yml b/.github/workflows/unfolding.yml index dbda27f8b..829a73bb2 100644 --- a/.github/workflows/unfolding.yml +++ b/.github/workflows/unfolding.yml @@ -264,7 +264,7 @@ jobs: scripts/ci/run_with_singularity.sh scripts/ci/setup_and_run_python.sh narf/scripts/fitting/combinetf2.py $WREMNANTS_OUTDIR/WMass_qGen_ptGen_absEtaGen_unfolding/WMass.hdf5 -o $WREMNANTS_OUTDIR/WMass_qGen_ptGen_absEtaGen_unfolding/fitresults.hdf5 --externalPostfit $WREMNANTS_OUTDIR/WMass_eta_pt_charge_unfolding/fitresults.hdf5 - --saveHists --computeHistErrors --project ch0 ptGen --project ch0 absEtaGen + --saveHists --computeHistErrors --computeHistCov --computeHistImpacts --noChi2 --project ch0 ptGen --project ch0 absEtaGen - name: wmass unfolding plot xsec run: >- @@ -298,7 +298,7 @@ jobs: - name: wmass theoryfit combine setup run: >- scripts/ci/run_with_singularity.sh scripts/ci/setup_and_run_python.sh scripts/combine/setupCombine.py -i $HIST_FILE - --fitresult ${WREMNANTS_OUTDIR}/WMass_eta_pt_charge_unfolding/fitresults.hdf5 --fitvar qGen-ptGen-absEtaGen + --fitresult ${WREMNANTS_OUTDIR}/WMass_qGen_ptGen_absEtaGen_unfolding/fitresults.hdf5 --fitvar qGen-ptGen-absEtaGen -o $WREMNANTS_OUTDIR --postfix theoryfit --baseName xnorm --realData - name: wmass combine fit @@ -357,7 +357,7 @@ jobs: - name: wlike analysis run: >- scripts/ci/run_with_singularity.sh scripts/ci/setup_and_run_python.sh scripts/histmakers/mz_wlike_with_mu_eta_pt.py --dataPath $DATAPATH - --analysisMode unfolding -j $((NTHREADS)) --maxFiles $((MAX_FILES)) --forceDefaultName -o $WREMNANTS_OUTDIR --postfix unfolding + --analysisMode unfolding --poiAsNoi -j $((NTHREADS)) --maxFiles $((MAX_FILES)) --forceDefaultName -o $WREMNANTS_OUTDIR --postfix unfolding wlike-unfolding: @@ -386,13 +386,13 @@ jobs: - name: wlike plot response matrix run: >- scripts/ci/run_with_singularity.sh scripts/ci/setup_and_run_python.sh scripts/plotting/response_matrix.py - --axes "pt-ptGen" "abs(eta)-absEtaGen" --procFilters Zmumu -p mz_wlike -o $WEB_DIR -f $PLOT_DIR $HIST_FILE + --axes "pt-ptGen" "abs(eta)-absEtaGen" --procFilters Zmumu -p mz_wlike -o $WEB_DIR -f $PLOT_DIR $HIST_FILE --histName yieldsUnfolding - name: wlike combine unfolding setup run: >- scripts/ci/run_with_singularity.sh scripts/ci/setup_and_run_python.sh scripts/combine/setupCombine.py - --analysisMode unfolding -i $HIST_FILE --lumiScale $LUMI_SCALE -o $WREMNANTS_OUTDIR --postfix unfolding - --genAxes qGen-ptGen-absEtaGen --sparse + --analysisMode unfolding --poiAsNoi -i $HIST_FILE --lumiScale $LUMI_SCALE -o $WREMNANTS_OUTDIR --postfix unfolding + --scaleNormXsecHistYields '0.01' --genAxes ptGen-absEtaGen-qGen - name: wlike combine unfolding fit run: >- @@ -403,21 +403,21 @@ jobs: - name: wlike combine xsec setup run: >- scripts/ci/run_with_singularity.sh scripts/ci/setup_and_run_python.sh scripts/combine/setupCombine.py - --analysisMode unfolding -i $HIST_FILE --lumiScale $LUMI_SCALE -o $WREMNANTS_OUTDIR --postfix unfolding + --analysisMode unfolding --poiAsNoi -i $HIST_FILE --lumiScale $LUMI_SCALE -o $WREMNANTS_OUTDIR --postfix unfolding --baseName xnorm --filterProcGroups Zmumu - --fitvar qGen-ptGen-absEtaGen --genAxes qGen-ptGen-absEtaGen + --scaleNormXsecHistYields '0.01' --fitvar ptGen-absEtaGen-qGen --genAxes ptGen-absEtaGen-qGen - name: wlike combine unfolded xsec run: >- scripts/ci/run_with_singularity.sh scripts/ci/setup_and_run_python.sh narf/scripts/fitting/combinetf2.py - $WREMNANTS_OUTDIR/ZMassWLike_qGen_ptGen_absEtaGen_unfolding/ZMassWLike.hdf5 -o $WREMNANTS_OUTDIR/ZMassWLike_qGen_ptGen_absEtaGen_unfolding/fitresults.hdf5 + $WREMNANTS_OUTDIR/ZMassWLike_ptGen_absEtaGen_qGen_unfolding/ZMassWLike.hdf5 -o $WREMNANTS_OUTDIR/ZMassWLike_ptGen_absEtaGen_qGen_unfolding/fitresults.hdf5 --externalPostfit $WREMNANTS_OUTDIR/ZMassWLike_eta_pt_charge_unfolding/fitresults.hdf5 - --saveHists --computeHistErrors --project ch0 ptGen --project ch0 absEtaGen + --saveHists --computeHistErrors --computeHistCov --computeHistImpacts --noChi2 --project ch0 ptGen --project ch0 absEtaGen - name: wlike unfolding plot xsec run: >- scripts/ci/run_with_singularity.sh scripts/ci/setup_and_run_python.sh scripts/plotting/postfitPlots.py - ${WREMNANTS_OUTDIR}/ZMassWLike_qGen_ptGen_absEtaGen_unfolding/fitresults.hdf5 -o $WEB_DIR -f $PLOT_DIR + ${WREMNANTS_OUTDIR}/ZMassWLike_ptGen_absEtaGen_qGen_unfolding/fitresults.hdf5 -o $WEB_DIR -f $PLOT_DIR --rrange 0.0 2.0 --unfoldedXsec --noUncertainty --noChisq --legCols 1 --noSciy --logoPos 0 --yscale 1.3 --project ch0 ptGen --project ch0 absEtaGen @@ -449,13 +449,13 @@ jobs: - name: wlike theoryfit combine setup run: >- scripts/ci/run_with_singularity.sh scripts/ci/setup_and_run_python.sh scripts/combine/setupCombine.py -i $HIST_FILE - --fitresult ${WREMNANTS_OUTDIR}/ZMassWLike_eta_pt_charge_unfolding/fitresults.hdf5 --fitvar qGen-ptGen-absEtaGen + --fitresult ${WREMNANTS_OUTDIR}/ZMassWLike_ptGen_absEtaGen_qGen_unfolding/fitresults.hdf5 --fitvar qGen-ptGen-absEtaGen -o $WREMNANTS_OUTDIR --postfix theoryfit --baseName xnorm --realData - name: wlike combine fit run: >- scripts/ci/run_with_singularity.sh scripts/ci/setup_and_run_python.sh narf/scripts/fitting/combinetf2.py - $WREMNANTS_OUTDIR/ZMassWLike_qGen_ptGen_absEtaGen_theoryfit/WMass.hdf5 -o $WREMNANTS_OUTDIR/ZMassWLike_qGen_ptGen_absEtaGen_theoryfit/fitresults.hdf5 + $WREMNANTS_OUTDIR/ZMassWLike_qGen_ptGen_absEtaGen_theoryfit/ZMassWLike.hdf5 -o $WREMNANTS_OUTDIR/ZMassWLike_qGen_ptGen_absEtaGen_theoryfit/fitresults.hdf5 --doImpacts --globalImpacts --saveHists --computeHistErrors --chisqFit --externalCovariance -t 0 - name: wlike theoryfit impacts @@ -566,7 +566,7 @@ jobs: scripts/ci/run_with_singularity.sh scripts/ci/setup_and_run_python.sh narf/scripts/fitting/combinetf2.py $WREMNANTS_OUTDIR/ZMassDilepton_ptVGen_absYVGen_unfolding/ZMassDilepton.hdf5 -o $WREMNANTS_OUTDIR/ZMassDilepton_ptVGen_absYVGen_unfolding/fitresults.hdf5 --externalPostfit $WREMNANTS_OUTDIR/ZMassDilepton_ptll_yll_unfolding/fitresults.hdf5 - --saveHists --computeHistErrors + --saveHists --computeHistErrors --computeHistCov --computeHistImpacts --noChi2 --project ch0 ptVGen --project ch0 absYVGen - name: dilepton unfolding plot pulls run: >- @@ -575,10 +575,12 @@ jobs: output --outFolder $WEB_DIR/$PLOT_DIR -o pulls_unfolding_ptll.html -n 50 --otherExtensions png pdf - name: dilepton unfolding plot xsec - run: >- - scripts/ci/run_with_singularity.sh scripts/ci/setup_and_run_python.sh scripts/plotting/unfolding_xsec.py $WREMNANTS_OUTDIR/ZMassDilepton_ptll_yll_unfolding/fitresults.hdf5 - --histfile $HIST_FILE --varNames scetlib_dyturboCorr --varLabels SCETLib $Omega_\nu=0.5$ --selectAxis vars --selectEntries 'omega_nu0.5' - -o $WEB_DIR -f $PLOT_DIR -v 4 --rrange 0.9 1.1 --logy -t 'utilities/styles/nuisance_translate.json' --grouping max + run: >- + scripts/ci/run_with_singularity.sh scripts/ci/setup_and_run_python.sh scripts/plotting/postfitPlots.py + $WREMNANTS_OUTDIR/ZMassDilepton_ptll_yll_unfolding/fitresults.hdf5 -o $WEB_DIR -f $PLOT_DIR + --rrange 0.9 1.1 --unfoldedXsec --noUncertainty --noChisq --legCols 1 --noSciy --logoPos 0 --yscale 1.3 + --project ch0 ptVGen --project ch0 absYVGen + copy-clean: runs-on: [self-hosted, linux, x64] diff --git a/narf b/narf index cf0f9e584..b8b66e6e0 160000 --- a/narf +++ b/narf @@ -1 +1 @@ -Subproject commit cf0f9e5841d6a21f06b720a85a3b2a6d1b7cf94a +Subproject commit b8b66e6e0caf6005e4b8c44bf380c40c9c0704ff diff --git a/scripts/plotting/postfitPlots.py b/scripts/plotting/postfitPlots.py index dabbc77f2..7c7b76c44 100644 --- a/scripts/plotting/postfitPlots.py +++ b/scripts/plotting/postfitPlots.py @@ -247,6 +247,7 @@ def make_plot( labels=None, hup=None, hdown=None, + h_data_stat=None, variation="", suffix="", chi2=None, @@ -324,6 +325,10 @@ def make_plot( hdown = [ hh.unrolledHist(h, binwnorm=binwnorm, obs=axes_names) for h in hdown ] + if h_data_stat is not None: + h_data_stat = hh.unrolledHist( + h_data_stat, binwnorm=binwnorm, obs=axes_names + ) if args.normToData: scale = h_data.values().sum() / h_inclusive.values().sum() @@ -414,6 +419,24 @@ def make_plot( zorder=2, flow="none", ) + if h_data_stat is not None: + var_stat = h_data_stat.values() ** 2 + h_data_stat = h_data.copy() + h_data_stat.variances()[...] = var_stat + + hep.histplot( + h_data_stat, + yerr=True, + histtype=histtype_data, + color="black", + # label=args.dataName, + binwnorm=binwnorm, + capsize=2, + ax=ax1, + alpha=1.0, + zorder=2, + flow="none", + ) if args.unfoldedXsec: hep.histplot( h_inclusive, @@ -483,6 +506,10 @@ def make_plot( if diff: h1 = hh.addHists(h_inclusive, h_inclusive, scale2=-1) h2 = hh.addHists(h_data, h_inclusive, scale2=-1) + if h_data_stat is not None: + h2_stat = hh.divideHists( + h_data_stat, h_inclusive, cutoff=0.01, rel_unc=True + ) else: h1 = hh.divideHists( h_inclusive, @@ -493,6 +520,10 @@ def make_plot( by_ax_name=False, ) h2 = hh.divideHists(h_data, h_inclusive, cutoff=0.01, rel_unc=True) + if h_data_stat is not None: + h2_stat = hh.divideHists( + h_data_stat, h_inclusive, cutoff=0.01, rel_unc=True + ) hep.histplot( h1, @@ -516,6 +547,17 @@ def make_plot( ax=ax2, flow="none", ) + if h_data_stat is not None: + hep.histplot( + h2_stat, + histtype="errorbar", + color="black", + yerr=True, + linewidth=2, + capsize=2, + ax=ax2, + flow="none", + ) # for uncertaity bands edges = h_inclusive.axes[0].edges @@ -718,7 +760,9 @@ def make_plots( labels, colors, hist_var=None, + hist_data_stat=None, channel="", + lumi=1, *opts, **kwopts, ): @@ -730,6 +774,14 @@ def make_plots( l if p not in args.suppressProcsLabel else None for l, p in zip(labels, procs) ] + if args.unfoldedXsec: + # convert number in cross section in pb + hist_data = hh.scaleHist(hist_data, 1.0 / (lumi * 1000)) + hist_inclusive = hh.scaleHist(hist_inclusive, 1.0 / (lumi * 1000)) + hist_stack = [hh.scaleHist(h, 1.0 / (lumi * 1000)) for h in hist_stack] + if hist_data_stat is not None: + hist_data_stat = hh.scaleHist(hist_data_stat, 1.0 / (lumi * 1000)) + if hist_var is not None: hists_down = [ hist_var[{"downUpVar": 0, "vars": n}].project(*[a.name for a in axes]) @@ -739,6 +791,8 @@ def make_plots( hist_var[{"downUpVar": 1, "vars": n}].project(*[a.name for a in axes]) for n in varNames ] + hists_up = [hh.scaleHist(h, 1.0 / (lumi * 1000)) for h in hists_up] + hists_down = [hh.scaleHist(h, 1.0 / (lumi * 1000)) for h in hists_down] else: hists_down = None hists_up = None @@ -765,6 +819,10 @@ def make_plots( h_data = hist_data[idxs] h_inclusive = hist_inclusive[idxs] h_stack = [h[idxs] for h in hist_stack] + if hist_data_stat is not None: + h_data_stat = hist_data_stat[idxs] + else: + h_data_stat = None if hist_var is not None: hdown = [h[idxs] for h in hists_down] @@ -806,7 +864,9 @@ def make_plots( suffix=suffix, hup=hup, hdown=hdown, + h_data_stat=h_data_stat, selection=selection, + lumi=lumi, *opts, **kwopts, ) @@ -821,6 +881,8 @@ def make_plots( suffix=channel, hup=hists_up, hdown=hists_down, + h_data_stat=hist_data_stat, + lumi=lumi, *opts, **kwopts, ) @@ -863,15 +925,20 @@ def make_plots( logger.info(f"Skip masked channel {channel}") continue - if not args.unfoldedXsec: + hist_data_stat = None + if args.unfoldedXsec: + hist_data = fitresult[f"hist_{fittype}_inclusive"][channel].get() + if f"hist_global_impacts_grouped_{fittype}_inclusive" in fitresult.keys(): + hist_data_stat = fitresult[ + f"hist_global_impacts_grouped_{fittype}_inclusive" + ][channel].get()[{"impacts": "stat"}] + hist_inclusive = fitresult[f"hist_prefit_inclusive"][channel].get() + hist_stack = [] + else: hist_data = fitresult["hist_data_obs"][channel].get() hist_inclusive = fitresult[f"hist_{fittype}_inclusive"][channel].get() hist_stack = fitresult[f"hist_{fittype}"][channel].get() hist_stack = [hist_stack[{"processes": p}] for p in procs] - else: - hist_data = fitresult[f"hist_prefit_inclusive"][channel].get() - hist_inclusive = fitresult[f"hist_{fittype}_inclusive"][channel].get() - hist_stack = [] # vary poi by postfit uncertainty if varNames is not None: @@ -915,19 +982,13 @@ def make_plots( h.values()[...] = 1e5 * np.log(h.values()) h.variances()[...] = 1e10 * (h.variances()) / np.square(or_vals) - lumi = info["lumi"] - if args.unfoldedXsec: - # convert number in cross section in pb - hist_data = hh.scaleHist(hist_data, 1.0 / (lumi * 1000)) - hist_inclusive = hh.scaleHist(hist_inclusive, 1.0 / (lumi * 1000)) - hist_stack = [hh.scaleHist(h, 1.0 / (lumi * 1000)) for h in hist_stack] - make_plots( hist_data, hist_inclusive, hist_stack, info["axes"], hist_var=hist_var, + hist_data_stat=hist_data_stat, channel=channel, procs=procs, labels=labels, @@ -935,17 +996,22 @@ def make_plots( chi2=chi2, saturated_chi2=saturated_chi2, meta=meta, - lumi=lumi, + lumi=info["lumi"], ) else: - for result in fitresult["projections"]: - channel = result["channel"] - axes = result["axes"] - - for projection in args.project: - if channel == projection[0] and set(axes) == set(projection[1:]): + for projection in args.project: + channel = projection[0] + axes = projection[1:] + + result = None + for fr in fitresult["projections"]: + if channel == fr["channel"] and set(axes) == set(fr["axes"]): + result = fr break else: + logger.warning( + f"Projection {projection} not found in fitresult, available projection are: {[(f['channel'], f['axes']) for f in fitresult['projections']]}. Skip!" + ) continue axes = [ @@ -957,15 +1023,20 @@ def make_plots( else: chi2 = None - if not args.unfoldedXsec: + hist_data_stat = None + if args.unfoldedXsec: + hist_data = result[f"hist_{fittype}_inclusive"].get() + if f"hist_global_impacts_grouped_{fittype}_inclusive" in result.keys(): + hist_data_stat = result[ + f"hist_global_impacts_grouped_{fittype}_inclusive" + ].get()[{"impacts": "stat"}] + hist_inclusive = result[f"hist_prefit_inclusive"].get() + hist_stack = [] + else: hist_data = result["hist_data_obs"].get() hist_inclusive = result[f"hist_{fittype}_inclusive"].get() hist_stack = result[f"hist_{fittype}"].get() hist_stack = [hist_stack[{"processes": p}] for p in procs] - else: - hist_data = result[f"hist_prefit_inclusive"].get() - hist_inclusive = result[f"hist_{fittype}_inclusive"].get() - hist_stack = [] # vary poi by postfit uncertainty if varNames is not None: @@ -973,27 +1044,20 @@ def make_plots( else: hist_var = None - lumi = meta_input["channel_info"][channel]["lumi"] - - if args.unfoldedXsec: - # convert number in cross section in pb - hist_data = hh.scaleHist(hist_data, 1.0 / (lumi * 1000)) - hist_inclusive = hh.scaleHist(hist_inclusive, 1.0 / (lumi * 1000)) - hist_stack = [hh.scaleHist(h, 1.0 / (lumi * 1000)) for h in hist_stack] - make_plots( hist_data, hist_inclusive, hist_stack, axes, hist_var=hist_var, + hist_data_stat=hist_data_stat, channel=channel, procs=procs, labels=labels, colors=colors, chi2=chi2, meta=meta, - lumi=lumi, + lumi=meta_input["channel_info"][channel]["lumi"], ) if output_tools.is_eosuser_path(args.outpath) and args.eoscp: