Merge pull request #415 from dbruschi/main

Updated veto selection to use only global muons, implementation of Veto Scale Factors and Uncertainties

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"]

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)

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

wremnants/include/muon_efficiencies_veto.h

@@ -0,0 +1,149 @@
+#ifndef WREMNANTS_MUON_EFFICIENCIES_VETO_H
+#define WREMNANTS_MUON_EFFICIENCIES_VETO_H
+
+#include <boost/histogram/axis.hpp>
+#include <array>
+#include "defines.h"
+
+namespace wrem {
+
+	template<typename HIST_SF, int NEtaBins, int NPtEigenBins, int NSysts, int Steps>
+	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<const HIST_SF>(std::move(sf_veto_plus))),
+			sf_veto_minus_(std::make_shared<const HIST_SF>(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<double, Eigen::Sizes<Steps, NSysts>>;
+
+		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<double, Eigen::Sizes<NEtaBins, NPtEigenBins, 2>>;
+
+		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<const HIST_SF> sf_veto_plus_;
+		std::shared_ptr<const HIST_SF> sf_veto_minus_;
+		double minpt_;
+		double maxpt_;
+		double mineta_;
+		double maxeta_;
+
+	};
+
+	template<typename HIST_SF, int NEtaBins, int NPtEigenBins, int NSysts, int Steps>
+	class muon_efficiency_veto_helper :
+		public muon_efficiency_veto_helper_base<HIST_SF, NEtaBins, NPtEigenBins, NSysts, Steps> {
+
+	public:
+
+		using base_t = muon_efficiency_veto_helper_base<HIST_SF, NEtaBins, NPtEigenBins, NSysts, Steps>;
+
+		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<typename HIST_SF, int NEtaBins, int NPtEigenBins, int NSysts, int Steps>
+	class muon_efficiency_veto_helper_syst :
+		public muon_efficiency_veto_helper_base<HIST_SF, NEtaBins, NPtEigenBins, NSysts, Steps> {
+
+	public:
+
+		using base_t = muon_efficiency_veto_helper_base<HIST_SF, NEtaBins, NPtEigenBins, NSysts, Steps>;
+		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<typename HIST_SF, int NEtaBins, int NPtEigenBins, int NSysts, int Steps>
+	class muon_efficiency_veto_helper_stat :
+		public muon_efficiency_veto_helper_base<HIST_SF, NEtaBins, NPtEigenBins, NSysts, Steps> {
+
+	public:
+
+		using base_t = muon_efficiency_veto_helper_base<HIST_SF, NEtaBins, NPtEigenBins, NSysts, Steps>;
+		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

wremnants/include/utils.h

@@ -90,6 +90,14 @@ double deltaR2(float eta1, float phi1, float eta2, float phi2) {
     return deta*deta + dphi*dphi;
 }
 
+RVec<double> vectDeltaR2(RVec<float> eta1, RVec<float> phi1, RVec<float> eta2, RVec<float> phi2) {
+	RVec<double> 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<Int_t> 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<float,float> i,std::pair<float,float> j) { return (i.first>j.first); }
+
+float unmatched_postfsrMuon_var(const Vec_f& var, const Vec_f& pt, int hasMatchDR2idx) {
+
+  std::vector<std::pair<float,float> > 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<float,int> i,std::pair<float,int> j) { return (i.first>j.first); }
+
+int unmatched_postfsrMuon_charge(const Vec_i& var, const Vec_f& pt, int hasMatchDR2idx) {
+
+  std::vector<std::pair<float,int> > 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<int> postFSRLeptonsIdx(RVec<bool> postFSRleptons) {
   RVec<int> v;
   for (unsigned int i=0; i<postFSRleptons.size(); i++) {