Skip to content

Commit

Permalink
Merge pull request #520 from cippy/validateIsoMtZ
Browse files Browse the repository at this point in the history
Update Wlike histmaker allowing iso-mt axes to be added to the list of nominal axes
  • Loading branch information
davidwalter2 authored Aug 10, 2024
2 parents f867dae + 4062efa commit 696350d
Show file tree
Hide file tree
Showing 2 changed files with 162 additions and 4 deletions.
136 changes: 136 additions & 0 deletions scripts/analysisTools/w_mass_13TeV/validateIsoMtRegionsWithZ.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,136 @@
#!/usr/bin/env python3

from wremnants.datasets.datagroups import Datagroups
from wremnants import histselections as sel
#from wremnants import plot_tools,theory_tools,syst_tools
from utilities import boostHistHelpers as hh
from utilities import common, logging
from utilities.io_tools import input_tools, output_tools

import narf
import wremnants
from wremnants import theory_tools,syst_tools
import hist

import numpy as np

import pickle
import lz4.frame

import argparse
import os
import shutil
import re

## safe batch mode
import sys
args = sys.argv[:]
sys.argv = ['-b']
import ROOT
sys.argv = args
ROOT.gROOT.SetBatch(True)
ROOT.PyConfig.IgnoreCommandLineOptions = True

from copy import *

from scripts.analysisTools.plotUtils.utility import *
from scripts.analysisTools.w_mass_13TeV.plotPrefitTemplatesWRemnants import plotPrefitHistograms

if __name__ == "__main__":
parser = common_plot_parser()
parser.add_argument("inputfile", type=str, nargs=1, help="Input file with histograms (pkl.lz4 or hdf5 file)")
parser.add_argument("outdir", type=str, nargs=1, help="Output folder")
parser.add_argument("-n", "--baseName", type=str, help="Histogram name in the file (it depends on what study you run)", default="nominal_testIsoMtFakeRegions")
parser.add_argument('-p','--processes', default=None, nargs='*', type=str,
help='Choose what processes to plot, otherwise all are done')
parser.add_argument("--isoBinEdges", type=int, nargs='+', default=[0, 0.15, 0.3], help="Bins to test (overflow added automatically)")
parser.add_argument("--mtBinEdges", type=float, nargs='+', default=[0,22,45], help="mT bins to test (overflow added automatically)")
parser.add_argument("--charge", type=int, default=1, choices=[-1, 0, 1], help="Triggering muon charge selection (0 integrates both charges)")
parser.add_argument("--rr", "--ratio-range", dest="ratioRange", default=(0.9,1.1), type=float, nargs=2, help="Range for ratio plot")
args = parser.parse_args()

logger = logging.setup_logger(os.path.basename(__file__), args.verbose)

charges = { -1 : "minus", 0: "both", 1 : "plus" }
chargeTag = charges[args.charge]
chargeBin = 0 if args.charge == -1 else 1
outdirTag = f"trigMuonCharge_{chargeTag}_nonTrigMuonPassIso/"
###############################################################################################
fname = args.inputfile[0]
outdir_original = f"{args.outdir[0]}/{outdirTag}/"
outdir = createPlotDirAndCopyPhp(outdir_original, eoscp=args.eoscp)

ROOT.TH1.SetDefaultSumw2()

canvas = ROOT.TCanvas("canvas", "", 800, 700)
cwide = ROOT.TCanvas("cwide","",2400,600)
adjustSettings_CMS_lumi()
canvas1D = ROOT.TCanvas("canvas1D", "", 800, 900)

groups = Datagroups(fname)
datasets = groups.getNames()
if args.processes is not None and len(args.processes):
datasets = list(filter(lambda x: x in args.processes, datasets))
else:
datasets = list(filter(lambda x: x != "QCD", datasets))

logger.debug(f"Using these processes: {datasets}")
inputHistName = args.baseName
groups.setNominalName(inputHistName)
groups.loadHistsForDatagroups(inputHistName, syst="", procsToRead=datasets, applySelection=False)
histInfo = groups.getDatagroups() # keys are same as returned by groups.getNames()
s = hist.tag.Slicer()
chargeSlice = s[chargeBin] if args.charge != 0 else s[::hist.sum]
eps = 0.0001
isoBinEdges = args.isoBinEdges
mtBinEdges = args.mtBinEdges

for jiso, iso in enumerate(isoBinEdges):
isoLow = round(iso,2)
isoLowStr = str(isoLow).replace(".","p")
maxJiso = len(isoBinEdges) - 1
isoHighStr = str(round(isoBinEdges[jiso+1],2)) if jiso < maxJiso else "Inf"
isoHighStr = isoHighStr.replace(".","p")
isoTag = f"iso{isoLowStr}To{isoHighStr}"
for jmt,mt in enumerate(mtBinEdges):
mtLowStr = str(round(mt))
maxJmt = len(mtBinEdges) - 1
mtHighStr = str(round(mtBinEdges[jmt+1])) if jmt < maxJmt else "Inf"
mtTag = f"mt{mtLowStr}To{mtHighStr}"
logger.info(f"Processing bin: {isoLowStr} < relIso < {isoHighStr} && {mtLowStr} < mT < {mtHighStr} GeV")
#
hdata2D = None
hmc2D = []
binTag = f"{isoTag}_{mtTag}"
outdirBin = f"{outdir}/{binTag}/"
createPlotDirAndCopyPhp(outdirBin, eoscp=args.eoscp)
for d in datasets:
logger.info(f"Running on process {d}")
hin = histInfo[d].hists[inputHistName]
# logger.debug(hin.axes)
###
# the edges we use might require integrating more bins, so can't just use the jxxx index
mtSlice = s[complex(0,mt+eps)::hist.sum]
if jmt < maxJmt:
mtSlice = s[complex(0,mt+eps):complex(0,mtBinEdges[jmt+1]+eps):hist.sum]
isoSlice = s[complex(0,iso+eps)::hist.sum]
if jiso < maxJiso:
isoSlice = s[complex(0,iso+eps):complex(0,isoBinEdges[jiso+1]+eps):hist.sum]
h = hin[{"mt" : mtSlice,
"relIso" : isoSlice,
"charge": chargeSlice}]
#logger.debug(h.axes)
if d =="Data":
hdata2D = narf.hist_to_root(h)
hdata2D.SetName(f"{d}_{chargeTag}")
hdata2D.SetTitle(f"{d} {binTag} {chargeTag}")
else:
hmc2D.append(narf.hist_to_root(h))
hmc2D[-1].SetName(f"{d}_{chargeTag}")
hmc2D[-1].SetTitle(f"{d} {binTag} {chargeTag}")
# end of process loop
plotPrefitHistograms(hdata2D, hmc2D, outdirBin, xAxisName="Triggering muon #eta", yAxisName="Triggering muon p_{T} (GeV)",
chargeLabel=chargeTag, canvas=canvas, canvasWide=cwide, canvas1D=canvas1D,
ratioRange=args.ratioRange, lumi=16.8)

copyOutputToEos(outdir, outdir_original, eoscp=args.eoscp)
30 changes: 26 additions & 4 deletions scripts/histmakers/mz_wlike_with_mu_eta_pt.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@

parser.add_argument("--mtCut", type=int, default=common.get_default_mtcut(analysis_label), help="Value for the transverse mass cut in the event selection") # 40 for Wmass, thus be 45 here (roughly half the boson mass)
parser.add_argument("--muonIsolation", type=int, nargs=2, default=[1,1], choices=[-1, 0, 1], help="Apply isolation cut to triggering and not-triggering muon (in this order): -1/1 for failing/passing isolation, 0 for skipping it")
parser.add_argument("--addIsoMtAxes", action="store_true", help="Add iso/mT axes to the nominal ones. It is for tests to get uncertainties (mainly from SF) versus iso-mT to validate the goodness of muon SF in the fake regions. Isolation (on triggering muon) and mT cut are automatically overridden.")
parser.add_argument("--validateVetoSF", action="store_true", help="Add histogram for validation of veto SF, loading all necessary helpers. This requires using the veto selection on the non-triggering muon, with reduced pt cut")
parser.add_argument("--useGlobalOrTrackerVeto", action="store_true", help="Use global-or-tracker veto definition and scale factors instead of global only")
parser.add_argument("--useRefinedVeto", action="store_true", help="Temporary option, it uses a different computation of the veto SF (only implemented for global muons)")
Expand All @@ -43,12 +44,16 @@

parser = common.set_parser_default(parser, "aggregateGroups", ["Diboson", "Top", "Wtaunu", "Wmunu"])
parser = common.set_parser_default(parser, "excludeProcs", ["QCD"])
if args.addIsoMtAxes:
parser = common.set_parser_default(parser, "muonIsolation", [0, 1])

if isTheoryAgnostic:
parser = common.set_parser_default(parser, "excludeFlow", True)
if args.genAbsYVbinEdges and any(x < 0.0 for x in args.genAbsYVbinEdges):
raise ValueError("Option --genAbsYVbinEdges requires all positive values. Please check")

args = parser.parse_args() # parse again or new defaults won't be propagated

if args.useRefinedVeto and args.useGlobalOrTrackerVeto:
raise NotImplementedError("Options --useGlobalOrTrackerVeto and --useRefinedVeto cannot be used together at the moment.")
if args.validateVetoSF:
Expand Down Expand Up @@ -98,6 +103,17 @@
nominal_axes = [axis_eta, axis_pt, common.axis_charge]
nominal_cols = ["trigMuons_eta0", "trigMuons_pt0", "trigMuons_charge0"]

# for isoMt region validation and related tests
# use very high upper edge as a proxy for infinity (cannot exploit overflow bins in the fit)
# can probably use numpy infinity, but this is compatible with the root conversion
axis_mtCat = hist.axis.Variable([0, int(args.mtCut/2.), args.mtCut, 1000], name = "mt", underflow=False, overflow=False)
axis_isoCat = hist.axis.Variable([0, 0.15, 0.3, 100], name = "relIso",underflow=False, overflow=False)

nominal_axes = [axis_eta, axis_pt, common.axis_charge]
nominal_cols = ["trigMuons_eta0", "trigMuons_pt0", "trigMuons_charge0"]
if args.addIsoMtAxes:
nominal_axes.extend([axis_mtCat, axis_isoCat])
nominal_cols.extend(["transverseMass", "trigMuons_relIso0"])

if isUnfolding:
template_wpt = (template_maxpt-template_minpt)/args.genBins[0]
Expand Down Expand Up @@ -265,8 +281,10 @@ def build_graph(df, dataset):
if not passIsoBoth:
df = muon_selections.apply_iso_muons(df, args.muonIsolation[0], args.muonIsolation[1], isoBranch, isoThreshold)

df = df.Define("trigMuons_passIso0", f"{isoBranch}[trigMuons][0] < {isoThreshold}")
df = df.Define("nonTrigMuons_passIso0", f"{isoBranch}[nonTrigMuons][0] < {isoThreshold}")
df = df.Define("trigMuons_relIso0", f"{isoBranch}[trigMuons][0]")
df = df.Define("nonTrigMuons_relIso0", f"{isoBranch}[nonTrigMuons][0]")
df = df.Define("trigMuons_passIso0", f"trigMuons_relIso0 < {isoThreshold}")
df = df.Define("nonTrigMuons_passIso0", f"nonTrigMuons_relIso0 < {isoThreshold}")

df = muon_selections.select_z_candidate(df, mass_min, mass_max)

Expand Down Expand Up @@ -364,14 +382,17 @@ def build_graph(df, dataset):

df = df.Define("passWlikeMT", f"transverseMass >= {mtw_min}")

if not args.onlyMainHistograms:
if not args.onlyMainHistograms and not isUnfolding and not args.addIsoMtAxes:
axis_mt_coarse = hist.axis.Variable([0.0, 15.0, 30.0, 45.0, 60.0, 75.0, 90.0, 120.0], name = "mt", underflow=False, overflow=True)
axis_trigPassIso = hist.axis.Boolean(name = f"trig_passIso")
axis_nonTrigPassIso = hist.axis.Boolean(name = f"nonTrig_passIso")

nominal_bin = df.HistoBoost("nominal_isoMtBins", [*axes, axis_trigPassIso, axis_nonTrigPassIso, axis_mt_coarse], [*cols, "trigMuons_passIso0", "nonTrigMuons_passIso0", "transverseMass", "nominal_weight"])
results.append(nominal_bin)

nominal_testIsoMtFakeRegions = df.HistoBoost("nominal_testIsoMtFakeRegions", [*axes, axis_isoCat, axis_mtCat], [*cols, "trigMuons_relIso0", "transverseMass", "nominal_weight"])
results.append(nominal_testIsoMtFakeRegions)

axis_eta_nonTrig = hist.axis.Regular(template_neta, template_mineta, template_maxeta, name = "etaNonTrig", overflow=False, underflow=False)
ptMin_nonTrig = args.vetoRecoPt if args.validateVetoSF else template_minpt
nPtBins_nonTrig = round(template_maxpt - args.vetoRecoPt) if args.validateVetoSF else template_npt
Expand All @@ -393,7 +414,8 @@ def build_graph(df, dataset):
df = df.Define("deltaPhiMuonMet", "std::abs(wrem::deltaPhi(trigMuons_phi0,met_wlike_TV2.Phi()))")
df = df.Filter(f"deltaPhiMuonMet > {args.dphiMuonMetCut*np.pi}")

df = df.Filter("passWlikeMT")
if not args.addIsoMtAxes:
df = df.Filter("passWlikeMT")

if not args.onlyMainHistograms:
# plot reco vertex distribution before and after PU reweigthing
Expand Down

0 comments on commit 696350d

Please sign in to comment.