diff --git a/config/analysis.yml b/config/analysis.yml index aea3f78b66d3fcfb15da95f9aae35a1d7ff5df7f..45d81aa2e15fad9037809beae0e992695538de90 100644 --- a/config/analysis.yml +++ b/config/analysis.yml @@ -40,6 +40,7 @@ plotIt: line-width: 1 line-color: 1 line-style: 1 + order: 0 ttB_OOA: legend: 't#bar{t}+B OOA' mpl_legend: '$\mathrm{t}\overline{\mathrm{t}}+B$ OOA' @@ -47,6 +48,7 @@ plotIt: line-width: 1 line-color: 1 line-style: 1 + order: 1 ttcc: legend: 't#bar{t}+C' mpl_legend: '$\mathrm{t}\overline{\mathrm{t}}+C$' @@ -54,6 +56,7 @@ plotIt: line-width: 1 line-color: 1 line-style: 1 + order: 2 ttjj: legend: 't#bar{t}+light' mpl_legend: '$\mathrm{t}\overline{\mathrm{t}}$+light' @@ -61,12 +64,14 @@ plotIt: line-width: 1 line-color: 1 line-style: 1 + order: 3 VJets: legend: 'V+jets' fill-color: "#EEEEEE" line-width: 1 line-color: 1 line-style: 1 + order: 8 ttX: legend: 't#bar{t}X' mpl_legend: '$\mathrm{t}\overline{\mathrm{t}}+X$' @@ -74,6 +79,7 @@ plotIt: line-width: 1 line-color: 1 line-style: 1 + order: 5 ttV: legend: 't#bar{t}V' mpl_legend: '$\mathrm{t}\overline{\mathrm{t}}\mathrm{V}$' @@ -81,6 +87,7 @@ plotIt: line-width: 1 line-color: 1 line-style: 1 + order: 7 ttH: legend: 't#bar{t}H' mpl_legend: '$\mathrm{t}\overline{\mathrm{t}}\mathrm{H}$' @@ -88,6 +95,7 @@ plotIt: line-width: 1 line-color: 1 line-style: 1 + order: 6 ST: legend: 'Single t' mpl_legend: 'Single $\mathrm{t}$' @@ -95,6 +103,7 @@ plotIt: line-width: 1 line-color: 1 line-style: 1 + order: 4 systematics: # this list is only used for control plots - muon_ID diff --git a/config/samples_template.yml b/config/samples_template.yml index e8da5dc6991b64c20cb52d6e1bbd0b5f47c214c2..f4ec6effbf7b2bb1bed973b006ac42e901d271d0 100644 --- a/config/samples_template.yml +++ b/config/samples_template.yml @@ -337,6 +337,7 @@ TTToHadronic_TuneCP5_13TeV-powheg-pythia8: signal_subprocesses: *tt_sig_subprocs signal_tag: "powheg_5FS" cross-section: &xs_tt_0l 377.96 + pdf_full: True generated-events: genEventSumw split: -30 @@ -381,6 +382,7 @@ TTbb_4f_TTToHadronic_TuneCP5-Powheg-Openloops-Pythia8: signal_tag: "powheg_4FS" cross-section: &xs_ttbb_0l 19.88 # 4FS #cross-section: 16.141 # for 5FS + pdf_full: True pdf_mc: True generated-events: genEventSumw split: -10 diff --git a/python/baseTtbbPlotter.py b/python/baseTtbbPlotter.py index 0b5e336aca9f920313bba7ca1757e6c460011083..4bf09ecb8b4563856467b2ca05bc1eca449db47a 100644 --- a/python/baseTtbbPlotter.py +++ b/python/baseTtbbPlotter.py @@ -43,6 +43,7 @@ class baseTtbbModule(NanoAODModule): parser.add_argument("--samples", nargs='*', required=True, help="Sample template YML file") parser.add_argument("--roc-corr", action="store_true", help="Enable muon Rochester correction") parser.add_argument("--pdf-mode", choices=["simple", "full"], default="simple", help="Full PDF uncertainties (100 histogram variations) or simple (event-based envelope) (only if systematics enabled) (default: %(default)s)") + parser.add_argument("--top-pt", action="store_true", help="Apply top pt reweighting (and add a 'noTopPt' variation where it is not applied)") def _loadSampleLists(self, analysisCfg): # fill sample template using JSON files @@ -126,7 +127,15 @@ class baseTtbbModule(NanoAODModule): noSel = noSel.refine("mcWeight", weight=tree.genWeight, autoSyst=sample_doSysts) if sample_doSysts: - noSel = utils.addTheorySystematics(self, sample, sampleCfg, tree, noSel, pdf_mode=self.args.pdf_mode) + noSel = utils.addTheorySystematics(self, sample, sampleCfg, tree, noSel, + pdf_mode=self.args.pdf_mode, + # only add additional top pt reweighting variation + # if we're not applying it everywhere + topPt=not self.args.top_pt) + + if "subprocess" in sampleCfg and self.args.top_pt: + # apply top pt reweighting for every variation + noSel = utils.applyTopPtReweighting(tree, noSel) if "subprocess" in sampleCfg and splitTT: logger.info(f"Adding ttbar category cuts for {sampleCfg['subprocess']}") diff --git a/python/controlPlotter.py b/python/controlPlotter.py index 90ee549d468ec8a1a9d3b2d5956e545efae1bd12..7e21f1fff7db940c1940363b988b9aca49c4ccea 100644 --- a/python/controlPlotter.py +++ b/python/controlPlotter.py @@ -78,10 +78,10 @@ class controlPlotter(recoBaseTtbbPlotter): plots += utils.makeMergedPlots( [(f"1mu_4j_2b", oneMu4Jet2BSel), (f"1ele_4j_2b", oneEle4Jet2BSel)], - f"1lep_4j_2b", "nBDeepFlavM", EqBin(5, 2, 7), title="Number of deepFlavourM b jets", var=op.rng_len(s.bJetsM)) + f"1lep_4j_2b", "nBDeepFlavM", EqBin(5, 2, 7), title="Number of medium b-tagged jets", var=op.rng_len(s.bJetsM)) plots += utils.makeMergedPlots( [(f"1mu_4j_2b", oneMu4Jet2BSel), (f"1ele_4j_2b", oneEle4Jet2BSel)], - f"1lep_4j_2b", "nBDeepFlavT", EqBin(6, 0, 6), title="Number of deepFlavourT b jets", var=op.rng_len(s.bJetsT)) + f"1lep_4j_2b", "nBDeepFlavT", EqBin(6, 0, 6), title="Number of tight b-tagged jets", var=op.rng_len(s.bJetsT)) plots += utils.makeMergedPlots( [(f"1mu_4j_2b", oneMu4Jet2BSel), (f"1ele_4j_2b", oneEle4Jet2BSel)], @@ -110,13 +110,13 @@ class controlPlotter(recoBaseTtbbPlotter): plots += utils.makeMergedPlots( [(f"1mu_5j_3b", oneMu5Jet3BSel), (f"1ele_5j_3b", oneEle5Jet3BSel)], - f"1lep_5j_3b", "nBDeepFlavM", EqBin(4, 3, 7), title="Number of deepFlavourM b jets", var=op.rng_len(s.bJetsM)) + f"1lep_5j_3b", "nBDeepFlavM", EqBin(4, 3, 7), title="Number of medium b-tagged jets", var=op.rng_len(s.bJetsM)) plots += utils.makeMergedPlots( [(f"1mu_5j_3b", oneMu5Jet3BSel), (f"1ele_5j_3b", oneEle5Jet3BSel)], - f"1lep_5j_3b", "nBDeepFlavT", EqBin(5, 0, 5), title="Number of deepFlavourT b jets", var=op.rng_len(s.bJetsT)) + f"1lep_5j_3b", "nBDeepFlavT", EqBin(5, 0, 5), title="Number of tight b-tagged jets", var=op.rng_len(s.bJetsT)) plots += utils.makeMergedPlots( [(f"1mu_5j_3b", oneMu5Jet3BSel), (f"1ele_5j_3b", oneEle5Jet3BSel)], - f"1lep_5j_3b", "nLightDeepFlavM", EqBin(6, 0, 6), title="Number of deepFlavourM light jets", var=op.rng_len(s.lightJetsM)) + f"1lep_5j_3b", "nLightDeepFlavM", EqBin(6, 0, 6), title="Number of light-tagged jets", var=op.rng_len(s.lightJetsM)) ##### Plots for ==1 lepton, >=6 jets, >= 2 b jets (ttjj level) ###### @@ -128,10 +128,10 @@ class controlPlotter(recoBaseTtbbPlotter): plots += utils.makeMergedPlots( [(f"1mu_6j_2b", oneMu6Jet2BSel), (f"1ele_6j_2b", oneEle6Jet2BSel)], - f"1lep_6j_2b", "nBDeepFlavM", EqBin(5, 2, 7), title="Number of deepFlavourM b jets", var=op.rng_len(s.bJetsM)) + f"1lep_6j_2b", "nBDeepFlavM", EqBin(5, 2, 7), title="Number of medium b-tagged jets", var=op.rng_len(s.bJetsM)) plots += utils.makeMergedPlots( [(f"1mu_6j_2b", oneMu6Jet2BSel), (f"1ele_6j_2b", oneEle6Jet2BSel)], - f"1lep_6j_2b", "nBDeepFlavT", EqBin(6, 0, 6), title="Number of deepFlavourT b jets", var=op.rng_len(s.bJetsT)) + f"1lep_6j_2b", "nBDeepFlavT", EqBin(6, 0, 6), title="Number of tight b-tagged jets", var=op.rng_len(s.bJetsT)) for sel,lep,name in [(oneMu6Jet2BSel, s.muon, f"1mu_6j_2b"), (oneEle6Jet2BSel, s.electron, f"1ele_6j_2b")]: plots += cp.makeLeptonPlots(sel, lep, name, binScaling=2) @@ -157,10 +157,10 @@ class controlPlotter(recoBaseTtbbPlotter): f"1lep_6j_4b", "nJets", EqBin(7, 6, 13), title="Number of jets", var=op.rng_len(s.cleanedJets)) plots += utils.makeMergedPlots( [(f"1mu_6j_4b", oneMu6Jet4BSel), (f"1ele_6j_4b", oneEle6Jet4BSel)], - f"1lep_6j_4b", "nBDeepFlavT", EqBin(6, 0, 6), title="Number of deepFlavourT b jets", var=op.rng_len(s.bJetsT)) + f"1lep_6j_4b", "nBDeepFlavT", EqBin(6, 0, 6), title="Number of tight b-tagged jets", var=op.rng_len(s.bJetsT)) plots += utils.makeMergedPlots( [(f"1mu_6j_4b", oneMu6Jet4BSel), (f"1ele_6j_4b", oneEle6Jet4BSel)], - f"1lep_6j_4b", "nLightDeepFlavM", EqBin(6, 0, 6), title="Number of deepFlavourM light jets", var=op.rng_len(s.lightJetsM)) + f"1lep_6j_4b", "nLightDeepFlavM", EqBin(6, 0, 6), title="Number of light-tagged jets", var=op.rng_len(s.lightJetsM)) plots += utils.makeMergedPlots( [(f"1mu_6j_4b", oneMu6Jet4BSel), (f"1ele_6j_4b", oneEle6Jet4BSel)], diff --git a/python/definitions.py b/python/definitions.py index 4db9007116a9c2974b8ff838d3d7a086e99abc5c..a75b615983c051dfc8fc705751329e5351b4b8bc 100644 --- a/python/definitions.py +++ b/python/definitions.py @@ -179,10 +179,10 @@ def eleDef(era): absEtaSC = op.abs(ele.eta + ele.deltaEtaSC) ele_pt = ele.pt # / ele.eCorr # uncomment to use uncalibrated electron pt if era == '2016ULpreVFP' or era == '2016ULpostVFP': - eraCut = op.AND(ele_pt > 29, op.abs(ele.eta) < 2.4) + eraCut = op.AND(ele_pt > 29, op.abs(ele.eta) < 2.5) if era == '2017UL' or era == '2018UL': eraCut = op.OR( - op.AND(ele_pt > 34., op.abs(ele.eta) < 2.4), + op.AND(ele_pt > 34., op.abs(ele.eta) < 2.5), op.AND(ele_pt > 30., op.abs(ele.eta) < 2.1), ) return op.AND( diff --git a/python/genDefinitions.py b/python/genDefinitions.py index 0181f84fe24da9c600b696a118a9d3874a7f071b..e39b8dd8c51ab0cf5930b4358065f3ed157ce187 100644 --- a/python/genDefinitions.py +++ b/python/genDefinitions.py @@ -8,28 +8,25 @@ import utils def muonDef(mu): return op.AND( mu.pt > 26., op.abs(mu.eta) < 2.4, - #ID - op.OR(mu.pdgId == 13, mu.pdgId == -13), + op.abs(mu.pdgId) == 13, ) def vetoMuonDef(mu): return op.AND( mu.pt > 15., op.abs(mu.eta) < 2.4, - #ID - op.OR(mu.pdgId == 13, mu.pdgId == -13), + op.abs(mu.pdgId) == 13, ) def eleDef(ele): return op.AND( - op.AND(ele.pt > 29, op.abs(ele.eta) < 2.4), - op.OR(ele.pdgId == 11, ele.pdgId == -11), + op.AND(ele.pt > 29, op.abs(ele.eta) < 2.5), + op.abs(ele.pdgId) == 11, ) def vetoEleDef(ele): return op.AND( ele.pt > 15., op.abs(ele.eta) < 2.5, - #ID - op.OR(ele.pdgId == 11, ele.pdgId == -11), + op.abs(ele.pdgId) == 11, ) def jetDef(jet): diff --git a/python/recoBaseTtbbPlotter.py b/python/recoBaseTtbbPlotter.py index d118ab82d6366761f76a49c55093de5257d92f2d..dc72c575bbedcba30ba711077a715c553f9368ab 100644 --- a/python/recoBaseTtbbPlotter.py +++ b/python/recoBaseTtbbPlotter.py @@ -116,9 +116,12 @@ class recoBaseTtbbPlotter(baseTtbbPlotter): "topmass": None, "type3": None, "jes": None, - "jer0": "jer", - "jer1": "jer", } + if self.args.no_split_jer: + systMapping["jer"] = "jer" + else: + systMapping["jer0"] = "jer" + systMapping["jer1"] = "jer" self.bTagWeightFun = defs.makeBtagSFWPs(era, noSel, wps=["M", "T"], decorr_wps=self.args.decorr_btag, full_scheme=True, full_scheme_mapping=systMapping, prepostVFPfix=self.args.btag_fix) diff --git a/python/unfoldingConfigWriter.py b/python/unfoldingConfigWriter.py index 5f2e50a6958d0a745a60cecdc337e7284d8b536c..3c85c827cf2acbf59d431557442dd1a797f1e4c7 100644 --- a/python/unfoldingConfigWriter.py +++ b/python/unfoldingConfigWriter.py @@ -81,6 +81,7 @@ class ConfigWriter(): "pretty_name" : obs.pretty_name(), "ROOT_pretty_name" : obs.ROOT_pretty_name(), "overflow" : obs.overflow(), + "logy" : obs.logy(), "phase_space" : self.get_phase_space(plot_names.gen_plot), "pretty_phase_space" : self.get_pretty_phase_space(plot_names.gen_plot) } diff --git a/python/unfoldingDefinitions.py b/python/unfoldingDefinitions.py index 18ba03b9bdc6a03ee45d1cdce57b76160f3600b9..1855c20dd243264bcf22653d8ea6f4383bf0980c 100644 --- a/python/unfoldingDefinitions.py +++ b/python/unfoldingDefinitions.py @@ -69,6 +69,9 @@ class Observable(): def overflow(self): return False + def logy(self): + return False + def plot_opts(self): return {} @@ -133,6 +136,10 @@ class Jet_Pt(Observable): def overflow(self): return True + def logy(self): + return True + + class Jet_AbsEta(Observable): def __init__(self, **kwargs): super().__init__(**kwargs) @@ -203,9 +210,6 @@ class Jet_dPhi(SingleObservable): def ROOT_pretty_name(self, i=None): return r'|#Delta#phi(lj^{extra}_{1},b_{soft})|' - def unit(self): - return 'GeV' - class Njet(SingleObservable): def __init__(self, nMin=5, nMax=10, **kwargs): super().__init__(**kwargs) @@ -402,6 +406,26 @@ class Average_DR(SingleObservable): return r'#Delta R_{bb}^{avg}' class Bjet_Pt(Jet_Pt): + def __init__(self, **kwargs): + super().__init__(**kwargs) + def _simple_name(self, i=None): + return f'bjet{self._idx(i)+1}_pt' + def pretty_name(self, i=None): + return r'$p_{\mathrm{T}}(\mathrm{b}_{' + str(self._idx(i)+1) + r'})$' + def ROOT_pretty_name(self, i=None): + return r'p_{T}(b_{' + str(self._idx(i)+1) + r'})' +class Bjet_Pt_5j_3b(Bjet_Pt): + def __init__(self, **kwargs): + super().__init__(**kwargs) + _reco_bins = [ + [], # 1st not measured + [], # 2nd not measured + [ 25, 40, 60, 80, 100, 125, 150, 175, 200 ] ] + _gen_bins = [[], [], + [ 25, 60, 100, 150, 200 ] ] + self.bin_dct['gen'] = [ VarBin( bin_edges) for bin_edges in _gen_bins ] + self.bin_dct['reco'] = [ VarBin( bin_edges) for bin_edges in _reco_bins ] +class Bjet_Pt_6j_4b(Bjet_Pt): def __init__(self, **kwargs): super().__init__(**kwargs) _reco_bins = [ @@ -415,14 +439,6 @@ class Bjet_Pt(Jet_Pt): self.bin_dct['gen'] = [ VarBin( bin_edges) for bin_edges in _gen_bins ] self.bin_dct['reco'] = [ VarBin( bin_edges) for bin_edges in _reco_bins ] - def _simple_name(self, i=None): - return f'bjet{self._idx(i)+1}_pt' - - def pretty_name(self, i=None): - return r'$p_{\mathrm{T}}(\mathrm{b}_{' + str(self._idx(i)+1) + r'})$' - - def ROOT_pretty_name(self, i=None): - return r'p_{T}(b_{' + str(self._idx(i)+1) + r'})' class Bjet_AbsEta(Jet_AbsEta): def __init__(self, **kwargs): @@ -449,22 +465,24 @@ class bJet_Ht(Jet_Ht): class bJet_Ht_5j_3b(bJet_Ht): def __init__(self, **kwargs): super().__init__(**kwargs) - _reco_bins = [ 150, 220, 270, 300, 340, 370, 410, 440, 480, 510, 580, 640, 700, 800, 1000 ] + _reco_bins = [ 90, 130, 150, 200, 250, 300, 350, 400, 450, 520, 600, 670, 750, 820, 900 ] _gen_bins = join_bins(_reco_bins) + _gen_bins[0] = 75 self.bin_dct['gen'] = [ VarBin(_gen_bins) ] self.bin_dct['reco'] = [ VarBin(_reco_bins) ] class bJet_Ht_6j_4b(bJet_Ht): def __init__(self, **kwargs): super().__init__(**kwargs) - _reco_bins = [ 200, 280, 330, 400, 450, 500, 550, 625, 700, 775, 850, 925, 1000, 1200, 1500 ] + _reco_bins = [ 120, 200, 250, 300, 350, 400, 450, 520, 600, 670, 750, 850, 1000] _gen_bins = join_bins(_reco_bins) + _gen_bins[0] = 100 self.bin_dct['gen'] = [ VarBin(_gen_bins) ] self.bin_dct['reco'] = [ VarBin(_reco_bins) ] class lightJet_Ht(Jet_Ht): def __init__(self, **kwargs): super().__init__(**kwargs) - _reco_bins = [ 60, 110, 150, 200, 250, 300, 350, 400, 450, 520, 600, 700, 800 ] + _reco_bins = [ 60, 120, 200, 275, 350, 425, 500, 575, 650, 725, 800 ] _gen_bins = join_bins(_reco_bins) self.bin_dct['gen'] = [ VarBin(_gen_bins) ] self.bin_dct['reco'] = [ VarBin(_reco_bins) ] diff --git a/python/unfoldingPlotter.py b/python/unfoldingPlotter.py index da51bab11376ec30ec051fc89b1bf580a1db1efd..d951a4e6a9fcac26c3dfbb1e628a912cb08552a8 100644 --- a/python/unfoldingPlotter.py +++ b/python/unfoldingPlotter.py @@ -215,20 +215,6 @@ class unfoldingPlotter(recoBaseTtbbPlotter, genBaseTtbbPlotter): era = sampleCfg["era"] plots = [] - # Define some observable groups (according to different kwargs they will receive, - # depending on the phase spaces they're used in) - njet_observables = [ unfDefs.Njet ] - nbjet_observables = [ unfDefs.NBTaggedJets ] - bjet_observables = [ unfDefs.Bjet_Pt, unfDefs.Bjet_AbsEta ] - bjetpair_observables = [ unfDefs.Average_DR ] - nonw_lightjet_observables = [ unfDefs.extraLightJet_Pt ] - dPhi_lightjet_observables = [ unfDefs.Jet_dPhi ] - lightjet_observables = [ unfDefs.lightJet_Ht ] - extra_jet_observables = [ unfDefs.Extra_Jet_Pt, unfDefs.Extra_Jet_AbsEta, - unfDefs.Extra_Jet_dR, unfDefs.Extra_Jet_Sum_Pt, - unfDefs.Extra_Jet_Sum_AbsEta, unfDefs.Extra_Jet_Mass ] - maxMbb_observables = [ unfDefs.Max_Mbb ] - jet_collections, \ bJet_collections, \ bJetPair_collections, \ @@ -273,34 +259,43 @@ class unfoldingPlotter(recoBaseTtbbPlotter, genBaseTtbbPlotter): # Light Jet Phase Spaces (6J3B/7J4B) if (gen_nJet == 6 and gen_nB == 3) or (gen_nJet == 7 and gen_nB == 4): - plots += unfDefs.addAllMeasurements(lightjet_observables, + plots += unfDefs.addAllMeasurements([unfDefs.lightJet_Ht], selection, lightJet_collections, obs_kwargs=obs_kwargs, aux=aux, cr_def=cr_def, config=self.configWriter) - plots += unfDefs.addAllMeasurements(nonw_lightjet_observables, + plots += unfDefs.addAllMeasurements([unfDefs.extraLightJet_Pt], selection, nonWLightJet_collections, obs_kwargs=obs_kwargs, aux=aux, cr_def=cr_def, config=self.configWriter) - plots += unfDefs.addAllMeasurements(dPhi_lightjet_observables, + plots += unfDefs.addAllMeasurements([unfDefs.Jet_dPhi], selection, nonWLightJet_collections, softB_collections, obs_kwargs=obs_kwargs, aux=aux, cr_def=cr_def, config=self.configWriter) # Primary Phase Spaces elif (gen_nJet == 5 and gen_nB == 3) or (gen_nJet == 6 and gen_nB == 4): - plots += unfDefs.addAllMeasurements(njet_observables, + plots += unfDefs.addAllMeasurements([unfDefs.Njet], selection, jet_collections, obs_kwargs=obs_kwargs | {'nMin': gen_nJet}, aux=aux, cr_def=cr_def, config=self.configWriter) - # b jet pt and eta only for 3rd b jet (5j3b) or 3rd and 4th b jets (6j4b) - plots += unfDefs.addAllMeasurements(bjet_observables, + # b jet eta only for 3rd b jet (5j3b) or 3rd and 4th b jets (6j4b) + plots += unfDefs.addAllMeasurements([unfDefs.Bjet_AbsEta], selection, bJet_collections, obs_kwargs=obs_kwargs | {"minIndex": 2, "maxIndex": gen_nB - 1}, aux=aux, cr_def=cr_def, config=self.configWriter) if gen_nJet == 5 and gen_nB == 3: - plots += unfDefs.addAllMeasurements([unfDefs.Jet_Ht, unfDefs.bJet_Ht_5j_3b], + plots += unfDefs.addAllMeasurements([unfDefs.Jet_Ht], selection, jet_collections, obs_kwargs=obs_kwargs, aux=aux, cr_def=cr_def, config=self.configWriter) + plots += unfDefs.addAllMeasurements([unfDefs.bJet_Ht_5j_3b], + selection, bJet_collections, + obs_kwargs=obs_kwargs, aux=aux, + cr_def=cr_def, config=self.configWriter) + # b jet pt only for 3rd b jet + plots += unfDefs.addAllMeasurements([unfDefs.Bjet_Pt_5j_3b], + selection, bJet_collections, + obs_kwargs=obs_kwargs | {"minIndex": 2, "maxIndex": 2}, aux=aux, + cr_def=cr_def, config=self.configWriter) # number of b jets, only in 5j3b - plots += unfDefs.addAllMeasurements(nbjet_observables, + plots += unfDefs.addAllMeasurements([unfDefs.NBTaggedJets], selection, jet_collections, obs_kwargs=obs_kwargs | { "working_point": "M", @@ -312,20 +307,34 @@ class unfoldingPlotter(recoBaseTtbbPlotter, genBaseTtbbPlotter): aux=aux, cr_def=cr_def, config=self.configWriter) if gen_nJet == 6 and gen_nB == 4: - plots += unfDefs.addAllMeasurements([unfDefs.Jet_Ht, unfDefs.bJet_Ht_6j_4b], + plots += unfDefs.addAllMeasurements([unfDefs.Jet_Ht], selection, jet_collections, obs_kwargs=obs_kwargs, aux=aux, cr_def=cr_def, config=self.configWriter) + plots += unfDefs.addAllMeasurements([unfDefs.bJet_Ht_6j_4b], + selection, bJet_collections, + obs_kwargs=obs_kwargs, aux=aux, + cr_def=cr_def, config=self.configWriter) + # b jet pt only for 3rd and 4th b jet + plots += unfDefs.addAllMeasurements([unfDefs.Bjet_Pt_6j_4b], + selection, bJet_collections, + obs_kwargs=obs_kwargs | {"minIndex": 2, "maxIndex": 3}, aux=aux, + cr_def=cr_def, config=self.configWriter) + + extra_jet_observables = [ unfDefs.Extra_Jet_Pt, unfDefs.Extra_Jet_AbsEta, + unfDefs.Extra_Jet_dR, unfDefs.Extra_Jet_Sum_Pt, + unfDefs.Extra_Jet_Sum_AbsEta, unfDefs.Extra_Jet_Mass ] + # extra jet pt and eta: for 1st and 2nd jet in the pair (duh) plots += unfDefs.addAllMeasurements(extra_jet_observables, selection, minDRbbPair_collections, obs_kwargs=obs_kwargs | {"maxIndex": 1}, aux=aux, cr_def=cr_def, config=self.configWriter) - plots += unfDefs.addAllMeasurements(maxMbb_observables, selection, + plots += unfDefs.addAllMeasurements([unfDefs.Max_Mbb], selection, maxMbbPair_collections, obs_kwargs=obs_kwargs, aux=aux, cr_def=cr_def, config=self.configWriter) - plots += unfDefs.addAllMeasurements(bjetpair_observables, selection, + plots += unfDefs.addAllMeasurements([unfDefs.Average_DR], selection, bJetPair_collections, obs_kwargs=obs_kwargs, aux=aux, cr_def=cr_def, config=self.configWriter) diff --git a/python/utils.py b/python/utils.py index fb6f6f40c49598449b0d666351651b51c60ad3ec..4827fc0f1561420e793a9ee81ac40a35bc89c1ec 100644 --- a/python/utils.py +++ b/python/utils.py @@ -149,36 +149,63 @@ def makeMergedPlots(categDef, newCat, name, binning, var=None, **kwargs): #### Common tasks (systematics, sample splittings...) +def getTopPtWeight(tree): + def top_pt_weight(pt): + return 0.103 * op.exp(-0.0118 * pt) - 0.000134 * pt + 0.973 + + lastCopy = op.select(tree.GenPart, lambda p: (p.statusFlags >> 13) & 1) + tops = op.select(lastCopy, lambda p: p.pdgId == 6) + antitops = op.select(lastCopy, lambda p: p.pdgId == -6) + weight = op.switch(op.AND(op.rng_len(tops) >= 1, op.rng_len(antitops) >= 1), + op.sqrt(top_pt_weight(tops[0].pt) * top_pt_weight(antitops[0].pt)), + 1.) + + return weight + + +def applyTopPtReweighting(tree, sel): + logger.info("Applying Top Pt reweighting everywhere, with a 'noTopPt' variation without it") + + return sel.refine("topPt", weight=op.systematic(getTopPtWeight(tree), noTopPt=op.c_float(1.))) + + import numpy as np -@ROOT.Numba.Declare(['RVec'], 'float') -def computeHessianPDFUncertainty(weights): +@ROOT.Numba.Declare(['RVec', 'float'], 'float') +def computeHessianPDFUncertainty(weights, SF=1.): if len(weights) < 2: return 0. weights = np.asarray(weights) - return np.sqrt(np.sum((weights[1:] - weights[0])**2)) + return np.sqrt(np.sum((SF * weights[1:] - weights[0])**2)) + NO_MESCALE_WGTS = ["ttbb_4fs_openloops_sherpa"] NO_PS_WGTS = ["tt_herwig7", "ttbb_4fs_openloops_sherpa"] -def addTheorySystematics(plotter, sample, sampleCfg, tree, noSel, qcdVar=True, PSISR=True, PSFSR=True, pdf_mode="simple"): + +def addTheorySystematics(plotter, sample, sampleCfg, tree, noSel, qcdVar=True, PSISR=True, PSFSR=True, pdf_mode="simple", topPt=True): signal_tag = sampleCfg.get("signal_tag") + SF = 1. # FxFx has a strange bug where ME and PDF (Hessian) weights are divided by 2 + if signal_tag == "ttjets_amcatnloFXFX": + SF = 2. + SF = op.c_float(SF) + ##### QCD muR/muF ###### if qcdVar: if hasattr(tree, "LHEScaleWeight") and signal_tag not in NO_MESCALE_WGTS: if plotter.qcdScaleVarMode == "separate": # for muF and muR separately logger.info("Adding separate muF/muR systematics") - noSel = noSel.refine("qcdMuF", weight=op.systematic(op.c_float(1.), name="qcdMuF", up=tree.LHEScaleWeight[5], down=tree.LHEScaleWeight[3])) - noSel = noSel.refine("qcdMuR", weight=op.systematic(op.c_float(1.), name="qcdMuR", up=tree.LHEScaleWeight[7], down=tree.LHEScaleWeight[1])) - noSel = noSel.refine("qcdMuRF", weight=op.systematic(op.c_float(1.), name="qcdMuRF", up=tree.LHEScaleWeight[8], down=tree.LHEScaleWeight[0])) + noSel = noSel.refine("qcdMuF", weight=op.systematic(op.c_float(1.), name="qcdMuF", up=SF * tree.LHEScaleWeight[5], down=SF * tree.LHEScaleWeight[3])) + noSel = noSel.refine("qcdMuR", weight=op.systematic(op.c_float(1.), name="qcdMuR", up=SF * tree.LHEScaleWeight[7], down=SF * tree.LHEScaleWeight[1])) + noSel = noSel.refine("qcdMuRF", weight=op.systematic(op.c_float(1.), name="qcdMuRF", up=SF * tree.LHEScaleWeight[8], down=SF * tree.LHEScaleWeight[0])) elif plotter.qcdScaleVarMode == "combined": # for taking envelope of 7-point variation qcdScaleVariations = plotter.getQcdScaleVariations(sample, sampleCfg) if qcdScaleVariations: logger.info("Adding 7-point muF/muR systematics") - qcdScaleVariations = { varNm: tree.LHEScaleWeight[varIdx] for (varIdx, varNm) in qcdScaleVariations } + qcdScaleVariations = { varNm: SF * tree.LHEScaleWeight[varIdx] for (varIdx, varNm) in qcdScaleVariations } qcdScaleSyst = op.systematic(op.c_float(1.), name="qcdScale", **qcdScaleVariations) noSel = noSel.refine("qcdScale", weight=qcdScaleSyst) else: @@ -212,14 +239,14 @@ def addTheorySystematics(plotter, sample, sampleCfg, tree, noSel, qcdVar=True, P if hasattr(tree, "LHEPdfWeight"): if pdf_mode == "full" and sampleCfg.get("pdf_full", False): logger.info("Adding full PDF systematics") - pdfVars = { f"pdf{i}": tree.LHEPdfWeight[i] for i in range(*nPdfVars) } - elif pdf_mode == "simple": + pdfVars = { f"pdf{i}": SF * tree.LHEPdfWeight[i] for i in range(*nPdfVars) } + else: logger.info("Adding simplified PDF systematics") if pdf_mc: pdfSigma = op.rng_stddev(tree.LHEPdfWeight) else: sigmaCalc = op.extMethod("Numba::computeHessianPDFUncertainty", returnType="float") - pdfSigma = sigmaCalc(tree.LHEPdfWeight) + pdfSigma = sigmaCalc(tree.LHEPdfWeight, SF) pdfVars = { "pdfup": tree.LHEPdfWeight[0] + pdfSigma, "pdfdown": tree.LHEPdfWeight[0] - pdfSigma } if pdf_mc: logger.info("This sample has MC PDF variations") @@ -229,8 +256,15 @@ def addTheorySystematics(plotter, sample, sampleCfg, tree, noSel, qcdVar=True, P else: logger.warning("LHEPdfWeight not present in tree, PDF systematics will not be added") + + ###### Top Pt ###### + if topPt and "subprocess" in sampleCfg: + logger.info("Adding top pt reweighting as additional variation") + noSel = noSel.refine("topPt", weight=op.systematic(op.c_float(1.), topPt=getTopPtWeight(tree))) + return noSel + def ttJetFlavCuts(subProc, tree): if subProc == "ttbb": return (tree.genTtbarId % 100) >= 53