diff --git a/narf b/narf index 542fc0f64..c2e5df7f7 160000 --- a/narf +++ b/narf @@ -1 +1 @@ -Subproject commit 542fc0f647f99731b6de712cb19dd8818a45c58e +Subproject commit c2e5df7f7ea2fcf0f638d6f01264ef5feb0658bb diff --git a/scripts/combine/setupCombine.py b/scripts/combine/setupCombine.py index b5b104d0f..29d47f2f7 100644 --- a/scripts/combine/setupCombine.py +++ b/scripts/combine/setupCombine.py @@ -686,6 +686,8 @@ def assertSample(name, startsWith=["W", "Z"], excludeMatch=[]): passToFakes=passSystToFakes) ## TODO: implement second lepton veto for low PU (both electrons and muons) if not lowPU: + pass + ''' # eta decorrelated nuisances decorrVarAxis = "eta" if "abseta" in fitvar: @@ -717,6 +719,7 @@ def assertSample(name, startsWith=["W", "Z"], excludeMatch=[]): passToFakes=passSystToFakes, scale=args.scaleZmuonVeto, ) + ''' else: cardTool.addLnNSystematic("CMS_background", processes=["Other"], size=1.15, group="CMS_background") @@ -774,6 +777,50 @@ def assertSample(name, startsWith=["W", "Z"], excludeMatch=[]): scale=scale, splitGroup=splitGroupDict, ) + + if wmass: + allEffTnP_veto = ["effStatTnP_veto_sf", "effSystTnP_veto"] + for name in allEffTnP_veto: + if "Syst" in name: + axes = ["veto_reco-veto_tracking-veto_idip", "n_syst_variations"] + axlabels = ["WPSYST", "_etaDecorr"] + nameReplace = [("WPSYST0", "reco"), ("WPSYST1", "tracking"), ("WPSYST2", "idip"), ("effSystTnP_veto", "effSyst_veto"), ("etaDecorr0", "fullyCorr") ] + scale = 1.0 + mirror = True + mirrorDownVarEqualToNomi=False + groupName = "muon_eff_veto_syst" + splitGroupDict = {f"{groupName}_{x}" : f".*effSyst_veto.*{x}" for x in list(["reco","tracking","idip"])} + else: + nameReplace = [] + mirror = True + mirrorDownVarEqualToNomi=False + if args.binnedScaleFactors: + axes = ["SF eta", "nPtBins", "SF charge"] + else: + axes = ["SF eta", "nPtEigenBins", "SF charge"] + axlabels = ["eta", "pt", "q"] + nameReplace = nameReplace + [("effStatTnP_veto_sf_", "effStat_veto_")] + scale = 1.0 + groupName = "muon_eff_veto_stat" + splitGroupDict = {} + if args.effStatLumiScale and "Syst" not in name: + scale /= math.sqrt(args.effStatLumiScale) + + cardTool.addSystematic( + name, + mirror=mirror, + mirrorDownVarEqualToNomi=mirrorDownVarEqualToNomi, + group=groupName, + systAxes=axes, + labelsByAxis=axlabels, + baseName=name+"_", + processes=['Zveto_samples'], + passToFakes=passSystToFakes, + systNameReplace=nameReplace, + scale=scale, + splitGroup=splitGroupDict, + ) + else: if datagroups.flavor in ["mu", "mumu"]: lepEffs = ["muSF_HLT_DATA_stat", "muSF_HLT_DATA_syst", "muSF_HLT_MC_stat", "muSF_HLT_MC_syst", "muSF_ISO_stat", "muSF_ISO_DATA_syst", "muSF_ISO_MC_syst", "muSF_IDIP_stat", "muSF_IDIP_DATA_syst", "muSF_IDIP_MC_syst"] diff --git a/scripts/histmakers/mw_with_mu_eta_pt.py b/scripts/histmakers/mw_with_mu_eta_pt.py index 8659d59e7..81bee159b 100644 --- a/scripts/histmakers/mw_with_mu_eta_pt.py +++ b/scripts/histmakers/mw_with_mu_eta_pt.py @@ -172,6 +172,7 @@ else: logger.info("Using smoothed scale factors and uncertainties") muon_efficiency_helper, muon_efficiency_helper_syst, muon_efficiency_helper_stat = wremnants.make_muon_efficiency_helpers_smooth(filename = args.sfFile, era = era, what_analysis = thisAnalysis, max_pt = axis_pt.edges[-1], isoEfficiencySmoothing = args.isoEfficiencySmoothing, smooth3D=args.smooth3dsf, isoDefinition=args.isolationDefinition) + muon_efficiency_veto_helper, muon_efficiency_veto_helper_syst, muon_efficiency_veto_helper_stat = wremnants.make_muon_efficiency_helpers_veto(era = era) logger.info(f"SF file: {args.sfFile}") @@ -400,6 +401,17 @@ def build_graph(df, dataset): df = df.Define("goodCleanJetsNoPt", "Jet_jetId >= 6 && (Jet_pt > 50 || Jet_puId >= 4) && abs(Jet_eta) < 2.4 && wrem::cleanJetsFromLeptons(Jet_eta,Jet_phi,Muon_correctedEta[vetoMuons],Muon_correctedPhi[vetoMuons],Electron_eta[vetoElectrons],Electron_phi[vetoElectrons])") df = df.Define("passIso", "goodMuons_relIso0 < 0.15") + ######################################################################## + # gen match to bare muons to select only prompt muons from MC processes, but also including tau decays (defined here because needed for veto SF) + if not dataset.is_data and not isQCDMC and not args.noGenMatchMC: + df = theory_tools.define_postfsr_vars(df) + df = df.Filter("wrem::hasMatchDR2(goodMuons_eta0,goodMuons_phi0,GenPart_eta[postfsrMuons],GenPart_phi[postfsrMuons],0.09)") + df = df.Define("postfsrMuons_inAcc", f"postfsrMuons && abs(GenPart_eta) < 2.4 && GenPart_pt > {args.vetoGenPartPt}") + df = df.Define("hasMatchDR2idx","wrem::hasMatchDR2idx(goodMuons_eta0,goodMuons_phi0,GenPart_eta[postfsrMuons_inAcc],GenPart_phi[postfsrMuons_inAcc],0.09)") + df = df.Define("unmatched_postfsrMuon_pt","wrem::unmatched_postfsrMuon_var(GenPart_pt[postfsrMuons_inAcc],GenPart_pt[postfsrMuons_inAcc],hasMatchDR2idx)") + df = df.Define("unmatched_postfsrMuon_eta","wrem::unmatched_postfsrMuon_var(GenPart_eta[postfsrMuons_inAcc],GenPart_pt[postfsrMuons_inAcc],hasMatchDR2idx)") + df = df.Define("GenPart_charge","-GenPart_pdgId[postfsrMuons_inAcc]/abs(GenPart_pdgId[postfsrMuons_inAcc])") + df = df.Define("unmatched_postfsrMuon_charge","wrem::unmatched_postfsrMuon_charge(GenPart_charge,GenPart_pt[postfsrMuons_inAcc],hasMatchDR2idx)") ######################################################################## # define event weights here since they are needed below for some helpers if dataset.is_data: @@ -427,7 +439,11 @@ def build_graph(df, dataset): if not isQCDMC and not args.noScaleFactors: df = df.Define("weight_fullMuonSF_withTrackingReco", muon_efficiency_helper, columnsForSF) weight_expr += "*weight_fullMuonSF_withTrackingReco" - + + if isZveto and not args.noGenMatchMC: + df = df.Define("weight_vetoSF_nominal", muon_efficiency_veto_helper, ["unmatched_postfsrMuon_pt","unmatched_postfsrMuon_eta","unmatched_postfsrMuon_charge"]) + weight_expr += "*weight_vetoSF_nominal" + logger.debug(f"Exp weight defined: {weight_expr}") df = df.Define("exp_weight", weight_expr) df = theory_tools.define_theory_weights_and_corrs(df, dataset.name, corr_helpers, args) @@ -436,7 +452,7 @@ def build_graph(df, dataset): df = theoryAgnostic_tools.define_helicity_weights(df) ######################################################################## - + ''' # gen match to bare muons to select only prompt muons from MC processes, but also including tau decays if not dataset.is_data and not isQCDMC and not args.noGenMatchMC: df = theory_tools.define_postfsr_vars(df) @@ -445,6 +461,7 @@ def build_graph(df, dataset): if isZveto: df = df.Define("postfsrMuons_inAcc", f"postfsrMuons && abs(GenPart_eta) < 2.4 && GenPart_pt > {args.vetoGenPartPt}") df = df.Define("ZvetoCondition", "Sum(postfsrMuons_inAcc) >= 2") + ''' if not args.noRecoil: leps_uncorr = ["Muon_pt[goodMuons][0]", "Muon_eta[goodMuons][0]", "Muon_phi[goodMuons][0]", "Muon_charge[goodMuons][0]"] @@ -609,11 +626,15 @@ def build_graph(df, dataset): if not isQCDMC and not args.noScaleFactors: df = syst_tools.add_muon_efficiency_unc_hists(results, df, muon_efficiency_helper_stat, muon_efficiency_helper_syst, axes, cols, what_analysis=thisAnalysis, smooth3D=args.smooth3dsf, storage_type=storage_type) + if isZveto and not args.noGenMatchMC: + df = syst_tools.add_muon_efficiency_veto_unc_hists(results, df, muon_efficiency_veto_helper_stat, muon_efficiency_veto_helper_syst, axes, cols, storage_type=storage_type) df = syst_tools.add_L1Prefire_unc_hists(results, df, muon_prefiring_helper_stat, muon_prefiring_helper_syst, axes, cols, storage_type=storage_type) # luminosity, as shape variation despite being a flat scaling to facilitate propagation to fakes df = syst_tools.add_luminosity_unc_hists(results, df, args, axes, cols, storage_type=storage_type) + ''' if isZveto: df = syst_tools.add_scaledByCondition_unc_hists(results, df, args, axes, cols, "weight_ZmuonVeto", "ZmuonVeto", "ZvetoCondition", 2.0, storage_type=storage_type) + ''' # n.b. this is the W analysis so mass weights shouldn't be propagated # on the Z samples (but can still use it for dummy muon scale) diff --git a/wremnants-data b/wremnants-data index ef5b5f397..619158845 160000 --- a/wremnants-data +++ b/wremnants-data @@ -1 +1 @@ -Subproject commit ef5b5f397a852d1144416b4a7f6c1a781c9b67e6 +Subproject commit 619158845203d60d447f2b5bae9e034a48d6c1f0 diff --git a/wremnants/__init__.py b/wremnants/__init__.py index cd2e021b7..f6e612ba3 100644 --- a/wremnants/__init__.py +++ b/wremnants/__init__.py @@ -15,6 +15,7 @@ from .muon_prefiring import make_muon_prefiring_helpers from .muon_efficiencies_smooth import make_muon_efficiency_helpers_smooth +from .muon_efficiencies_veto import make_muon_efficiency_helpers_veto from .muon_efficiencies_binned import make_muon_efficiency_helpers_binned from .muon_efficiencies_binned_vqt import make_muon_efficiency_helpers_binned_vqt from .muon_efficiencies_binned_vqt_integrated import make_muon_efficiency_helpers_binned_vqt_integrated diff --git a/wremnants/include/muon_efficiencies_veto.h b/wremnants/include/muon_efficiencies_veto.h new file mode 100644 index 000000000..d66b4685c --- /dev/null +++ b/wremnants/include/muon_efficiencies_veto.h @@ -0,0 +1,149 @@ +#ifndef WREMNANTS_MUON_EFFICIENCIES_VETO_H +#define WREMNANTS_MUON_EFFICIENCIES_VETO_H + +#include +#include +#include "defines.h" + +namespace wrem { + + template + class muon_efficiency_veto_helper_base { + + public: + muon_efficiency_veto_helper_base(HIST_SF &&sf_veto_plus, HIST_SF &&sf_veto_minus, double minpt, double maxpt, double mineta, double maxeta) : + sf_veto_plus_(std::make_shared(std::move(sf_veto_plus))), + sf_veto_minus_(std::make_shared(std::move(sf_veto_minus))), + minpt_(minpt), + maxpt_(maxpt), + mineta_(mineta), + maxeta_(maxeta) { + } + + double scale_factor(float pt, float eta, int charge) const { + + double sf = 1.0; + + if ((pt > minpt_) && (pt < maxpt_) && (eta > mineta_) && (eta < maxeta_)) { + auto const ix = sf_veto_plus_->template axis<0>().index(eta); + auto const iy = sf_veto_plus_->template axis<1>().index(pt); + if (charge>0) sf = sf_veto_plus_->at(ix,iy,0).value(); + else sf = sf_veto_minus_->at(ix,iy,0).value(); + } + + return sf; + + } + + using syst_tensor_t = Eigen::TensorFixedSize>; + + syst_tensor_t sf_syst_var(float pt, float eta, int charge) const { + + syst_tensor_t res; + res.setConstant(1.0); + + if ((pt > minpt_) && (pt < maxpt_) && (eta > mineta_) && (eta < maxeta_)) { + auto const ix = sf_veto_plus_->template axis<0>().index(eta); + auto const iy = sf_veto_plus_->template axis<1>().index(pt); + for (int step = 0; step < Steps; step++) { + for (int ns = 0; ns < NSysts; ns++) { + if ((ns == 0) || (ns == (ix+1))) { + if (charge>0) res(step, ns) = sf_veto_plus_->at(ix,iy,1+2*NPtEigenBins+step).value()/sf_veto_plus_->at(ix,iy,0).value(); + else res(step, ns) = sf_veto_minus_->at(ix,iy,1+2*NPtEigenBins+step).value()/sf_veto_minus_->at(ix,iy,0).value(); + } + else res(step, ns) = 1.; + } + } + } + return res; + + } + + using stat_tensor_t = Eigen::TensorFixedSize>; + + stat_tensor_t sf_stat_var(float pt, float eta, int charge) const { + stat_tensor_t res; + res.setConstant(1.0); + + if ((pt > minpt_) && (pt < maxpt_) && (eta > mineta_) && (eta < maxeta_)) { + auto const ix = sf_veto_plus_->template axis<0>().index(eta); + auto const iy = sf_veto_plus_->template axis<1>().index(pt); + for (int tensor_eigen_idx = 1; tensor_eigen_idx <= NPtEigenBins; tensor_eigen_idx++) { + if (charge > 0) res(ix, tensor_eigen_idx-1, 1) *= sf_veto_plus_->at(ix,iy,tensor_eigen_idx).value()/sf_veto_plus_->at(ix,iy,0).value(); + else res(ix, tensor_eigen_idx-1, 0) *= sf_veto_minus_->at(ix,iy,tensor_eigen_idx).value()/sf_veto_minus_->at(ix,iy,0).value(); + } + } + + return res; + } + + protected: + + std::shared_ptr sf_veto_plus_; + std::shared_ptr sf_veto_minus_; + double minpt_; + double maxpt_; + double mineta_; + double maxeta_; + + }; + + template + class muon_efficiency_veto_helper : + public muon_efficiency_veto_helper_base { + + public: + + using base_t = muon_efficiency_veto_helper_base; + + using base_t::base_t; + + muon_efficiency_veto_helper(const base_t &other) : base_t(other) {} + + double operator() (float pt, float eta, int charge) { + return base_t::scale_factor(pt, eta, charge); + } + + }; + + template + class muon_efficiency_veto_helper_syst : + public muon_efficiency_veto_helper_base { + + public: + + using base_t = muon_efficiency_veto_helper_base; + using tensor_t = typename base_t::syst_tensor_t; + + using base_t::base_t; + + muon_efficiency_veto_helper_syst(const base_t &other) : base_t(other) {} + + tensor_t operator() (float pt, float eta, int charge, double nominal_weight = 1.0) { + return nominal_weight*base_t::sf_syst_var(pt, eta, charge); + } + + }; + + template + class muon_efficiency_veto_helper_stat : + public muon_efficiency_veto_helper_base { + + public: + + using base_t = muon_efficiency_veto_helper_base; + using tensor_t = typename base_t::stat_tensor_t; + + using base_t::base_t; + + muon_efficiency_veto_helper_stat(const base_t &other) : base_t(other) {} + + tensor_t operator() (float pt, float eta, int charge, double nominal_weight = 1.0) { + return nominal_weight*base_t::sf_stat_var(pt, eta, charge); + } + + }; + +} + +#endif diff --git a/wremnants/include/utils.h b/wremnants/include/utils.h index f27909eb4..809a78a7c 100644 --- a/wremnants/include/utils.h +++ b/wremnants/include/utils.h @@ -90,6 +90,14 @@ double deltaR2(float eta1, float phi1, float eta2, float phi2) { return deta*deta + dphi*dphi; } +RVec vectDeltaR2(RVec eta1, RVec phi1, RVec eta2, RVec phi2) { + RVec vect; + for (unsigned int i=0U; i!=eta1.size(); i++) { + vect.push_back(deltaR2(eta1[i],phi1[i],eta2[i],phi2[i])); + } + return vect; +} + Vec_i cleanJetsFromLeptons(const Vec_f& Jet_eta, const Vec_f& Jet_phi, const Vec_f& Muon_eta, const Vec_f& Muon_phi, const Vec_f& Electron_eta, const Vec_f& Electron_phi) { Vec_i res(Jet_eta.size(), 1); // initialize to true and set to false whenever the jet overlaps with a muon @@ -217,6 +225,43 @@ RVec hasMatchDR2collWithSingle(const Vec_f &coll1_eta, const Vec_f &coll1 return resDR; } +int hasMatchDR2idx(const float& eta, const float& phi, const Vec_f& vec_eta, const Vec_f& vec_phi, const float dr2 = 0.09) { + + for (unsigned int jvec = 0; jvec < vec_eta.size(); ++jvec) { + if (deltaR2(eta, phi, vec_eta[jvec], vec_phi[jvec]) < dr2) return jvec; + } + return -1; + +} + +bool veto_vector_sorting(std::pair i,std::pair j) { return (i.first>j.first); } + +float unmatched_postfsrMuon_var(const Vec_f& var, const Vec_f& pt, int hasMatchDR2idx) { + + std::vector > ptsortingvector; + for (unsigned int i = 0; i < var.size(); i++) { + if (i!=hasMatchDR2idx) ptsortingvector.push_back(std::pair(pt[i],var[i])); + } + if (ptsortingvector.size() == 0) return -99; + std::sort(ptsortingvector.begin(),ptsortingvector.end(),veto_vector_sorting); + return ptsortingvector[0].second; + +} + +bool veto_vector_sorting_charge(std::pair i,std::pair j) { return (i.first>j.first); } + +int unmatched_postfsrMuon_charge(const Vec_i& var, const Vec_f& pt, int hasMatchDR2idx) { + + std::vector > ptsortingvector; + for (unsigned int i = 0; i < var.size(); i++) { + if (i!=hasMatchDR2idx) ptsortingvector.push_back(std::pair(pt[i],var[i])); + } + if (ptsortingvector.size() == 0) return -99; + std::sort(ptsortingvector.begin(),ptsortingvector.end(),veto_vector_sorting_charge); + return ptsortingvector[0].second; + +} + RVec postFSRLeptonsIdx(RVec postFSRleptons) { RVec v; for (unsigned int i=0; i {dlenSig} && SV_ntracks >= {ntracks}") return df +''' def select_veto_muons(df, nMuons=1, condition="==", ptCut=10.0, etaCut=2.4): # n.b. charge = -99 is a placeholder for invalid track refit/corrections (mostly just from tracks below @@ -24,6 +26,19 @@ def select_veto_muons(df, nMuons=1, condition="==", ptCut=10.0, etaCut=2.4): return df +''' + +def select_veto_muons(df, nMuons=1, condition="==", ptCut=15.0, etaCut=2.4): + + # n.b. charge = -99 is a placeholder for invalid track refit/corrections (mostly just from tracks below + # the pt threshold of 8 GeV in the nano production) + df = df.Define("vetoMuonsPre", "Muon_looseId && abs(Muon_dxybs) < 0.05 && Muon_correctedCharge != -99") + df = df.Define("vetoMuonsPre2", "vetoMuonsPre && Muon_isGlobal && Muon_highPurity && Muon_standalonePt > 15 && Muon_standaloneNumberOfValidHits > 0 && wrem::vectDeltaR2(Muon_standaloneEta, Muon_standalonePhi, Muon_correctedEta, Muon_correctedPhi) < 0.09") + df = df.Define("vetoMuons", f"vetoMuonsPre2 && Muon_correctedPt > {ptCut} && abs(Muon_correctedEta) < {etaCut}") + df = df.Filter(f"Sum(vetoMuons) {condition} {nMuons}") + + return df + def select_good_muons(df, ptLow, ptHigh, datasetGroup, nMuons=1, use_trackerMuons=False, use_isolation=False, isoDefinition="iso04", isoThreshold=0.15, condition="==", nonPromptFromSV=False, nonPromptFromLighMesonDecay=False): if use_trackerMuons: diff --git a/wremnants/syst_tools.py b/wremnants/syst_tools.py index c04b25ec6..a98616607 100644 --- a/wremnants/syst_tools.py +++ b/wremnants/syst_tools.py @@ -595,6 +595,23 @@ def add_muon_efficiency_unc_hists(results, df, helper_stat, helper_syst, axes, c return df +def add_muon_efficiency_veto_unc_hists(results, df, helper_stat, helper_syst, axes, cols, base_name="nominal", storage_type=hist.storage.Double()): + # TODO: update for dilepton + muon_columns_stat = ["unmatched_postfsrMuon_pt","unmatched_postfsrMuon_eta","unmatched_postfsrMuon_charge"] + muon_columns_syst = ["unmatched_postfsrMuon_pt","unmatched_postfsrMuon_eta","unmatched_postfsrMuon_charge"] + + df = df.Define("effStatTnP_veto_tensor", helper_stat, [*muon_columns_stat, "nominal_weight"]) + name = Datagroups.histName(base_name, syst="effStatTnP_veto_sf") + effStatTnP = df.HistoBoost(name, axes, [*cols, "effStatTnP_veto_tensor"], tensor_axes = helper_stat.tensor_axes, storage=storage_type) + results.append(effStatTnP) + + df = df.Define("effSystTnP_veto_weight", helper_syst, [*muon_columns_syst, "nominal_weight"]) + name = Datagroups.histName(base_name, syst="effSystTnP_veto") + effSystTnP = df.HistoBoost(name, axes, [*cols, "effSystTnP_veto_weight"], tensor_axes = helper_syst.tensor_axes, storage=storage_type) + results.append(effSystTnP) + + return df + def add_L1Prefire_unc_hists(results, df, helper_stat, helper_syst, axes, cols, base_name="nominal", addhelicity=False, storage_type=hist.storage.Double()): df = df.Define("muonL1PrefireStat_tensor", helper_stat, ["Muon_correctedEta", "Muon_correctedPt", "Muon_correctedPhi", "Muon_correctedCharge", "Muon_looseId", "nominal_weight"]) name = Datagroups.histName(base_name, syst=f"muonL1PrefireStat")