Skip to content

Commit 65db08f

Browse files
kyj519kyj519
authored
Merge pull request #52 from choij1589/main
Following PR from #50 - JES tag updates + etc.
2 parents 9698dff + 36b48fc commit 65db08f

32 files changed

Lines changed: 9144 additions & 201 deletions

‎AnalyzerTools/include/MyCorrection.h‎

Lines changed: 68 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
#include <vector>
77
#include <string>
88
#include <unordered_set>
9+
#include <variant>
910
using namespace std;
1011

1112
#include "TString.h"
@@ -51,7 +52,6 @@ class MyCorrection {
5152
inline float GetMuonISOSF(const TString &Muon_ISO_SF_Key, const Muon &muon, const variation syst = variation::nom, const TString &source = "") { return GetMuonIDSF(Muon_ISO_SF_Key, muon, syst); }
5253
inline float GetMuonTriggerSF(const TString &Muon_Trigger_SF_Key, const Muon &muon, const variation syst = variation::nom, const TString &source = "") { return GetMuonIDSF(Muon_Trigger_SF_Key, muon, syst); };
5354
inline float GetMuonISOSF(const TString &Muon_ISO_SF_Key, const RVec<Muon> &muons, const variation syst = variation::nom, const TString &source = "") { return GetMuonIDSF(Muon_ISO_SF_Key, muons, syst); }
54-
// float GetMuonTriggerSF(const TString &Muon_Trigger_SF_Key, const RVec<Muon> &muons, const variation syst = variation::nom, const TString &source = "");
5555
float GetMuonIDSF(const TString &Muon_ID_SF_Key, const Muon &muon, const variation syst = variation::nom) const;
5656
float GetMuonIDSF(const TString &Muon_ID_SF_Key, const RVec<Muon> &muons, const variation syst = variation::nom) const;
5757

@@ -129,6 +129,52 @@ class MyCorrection {
129129
// reweighting
130130
float GetTopPtReweight(const RVec<Gen> &gens) const;
131131

132+
// Safe evaluation function for correction sets with comprehensive error handling
133+
template<typename... Args>
134+
inline float safeEvaluate(const correction::Correction::Ref &cset,
135+
const string &function_name,
136+
const vector<correction::Variable::Type> &args) const {
137+
if (!cset) {
138+
cerr << "[MyCorrection::" << function_name << "] Error: Correction set is null" << endl;
139+
exit(EXIT_FAILURE);
140+
}
141+
142+
try {
143+
return cset->evaluate(args);
144+
} catch (const std::exception &e) {
145+
cerr << "[MyCorrection::" << function_name << "] Error during evaluation: " << e.what() << endl;
146+
cerr << "[MyCorrection::" << function_name << "] Arguments (" << args.size() << "): ";
147+
for (const auto &arg : args) {
148+
std::visit([](const auto &value) { cerr << value << " "; }, arg);
149+
}
150+
cerr << endl;
151+
exit(EXIT_FAILURE);
152+
}
153+
}
154+
155+
// Overload for CompoundCorrection
156+
template<typename... Args>
157+
inline float safeEvaluate(const correction::CompoundCorrection::Ref &cset,
158+
const string &function_name,
159+
const vector<correction::Variable::Type> &args) const {
160+
if (!cset) {
161+
cerr << "[MyCorrection::" << function_name << "] Error: CompoundCorrection set is null" << endl;
162+
exit(EXIT_FAILURE);
163+
}
164+
165+
try {
166+
return cset->evaluate(args);
167+
} catch (const std::exception &e) {
168+
cerr << "[MyCorrection::" << function_name << "] Error during evaluation: " << e.what() << endl;
169+
cerr << "[MyCorrection::" << function_name << "] Arguments (" << args.size() << "): ";
170+
for (const auto &arg : args) {
171+
std::visit([](const auto &value) { cerr << value << " "; }, arg);
172+
}
173+
cerr << endl;
174+
exit(EXIT_FAILURE);
175+
}
176+
}
177+
132178
private:
133179
struct EraConfig {
134180
string json_muon;
@@ -354,7 +400,27 @@ class MyCorrection {
354400
};
355401
return sys_string;
356402
};
357-
403+
404+
inline string getSystString_EGMScale(const variation syst) const {
405+
// Only for Run2
406+
if (! (Run == 2)) {
407+
throw runtime_error("[MyCorrection::getSystString_EGMScale] Use getSystString_EGM for Run3");
408+
}
409+
string sys_string = "";
410+
switch(syst) {
411+
case variation::nom:
412+
sys_string = "";
413+
break;
414+
case variation::up:
415+
sys_string = "scaleup";
416+
break;
417+
case variation::down:
418+
sys_string = "scaledown";
419+
break;
420+
};
421+
return sys_string;
422+
}
423+
358424
inline string getSystString_JME(const variation syst) const {
359425
string sys_string = "nom";
360426
switch (syst) {

‎AnalyzerTools/src/MyCorrection.cc‎

Lines changed: 119 additions & 134 deletions
Large diffs are not rendered by default.

‎Analyzers/include/DiLepton.h‎

Lines changed: 16 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -63,7 +63,7 @@ class DiLepton: public DiLeptonBase {
6363
void processWeightOnlySystematics(const Channel& channel, const Event& event, const RecoObjects& recoObjects, const RVec<Gen>& genParts);
6464

6565
// private methods
66-
Channel selectEvent(Event& ev, const RecoObjects& recoObjects);
66+
Channel selectEvent(Event& ev, const RecoObjects& recoObjects, const TString& syst);
6767
RecoObjects defineObjects(Event& ev, const RVec<Muon>& rawMuons,
6868
const RVec<Electron>& rawElectrons,
6969
const RVec<Jet>& rawJets,
@@ -80,6 +80,21 @@ class DiLepton: public DiLeptonBase {
8080
const RecoObjects& recoObjects,
8181
const WeightInfo& weights,
8282
const TString& syst = "Central");
83+
84+
// Cutflow functionality
85+
enum class CutStage {
86+
Initial = 0,
87+
NoiseFilter = 1,
88+
VetoMap = 2,
89+
LeptonSelection = 3,
90+
Trigger = 4,
91+
KinematicCuts = 5,
92+
JetRequirements = 6,
93+
BjetRequirements = 7,
94+
Final = 8
95+
};
96+
97+
void fillCutflow(CutStage stage, const Channel& channel, float weight, const TString& syst);
8398
};
8499

85100
#endif

‎Analyzers/src/AnalyzerCore.cc‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -837,7 +837,6 @@ RVec<Jet> AnalyzerCore::GetAllJets() {
837837
jet.SetPtEtaPhiM(correctedPt, Jet_eta[i], Jet_phi[i], correctedMass);
838838
jet.SetRawPt(rawPt);
839839
jet.SetOriginalPt(Jet_pt[i]);
840-
//jet.SetJESUncertainty(scaleUnc);
841840
jet.SetArea(Jet_area[i]);
842841
jet.SetOriginalIndex(i);
843842
jet.SetEnergyFractions(Jet_chHEF[i], Jet_neHEF[i], Jet_neEmEF[i], Jet_chEmEF[i], Jet_muEF[i]);

‎Analyzers/src/DiLepton.cc‎

Lines changed: 41 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -19,20 +19,22 @@ void DiLepton::initializeAnalyzer() {
1919
} else {
2020
systHelper = std::make_unique<SystematicHelper>(SKNANO_HOME + "/AnalyzerTools/DiLeptonSystematic.yaml", MCSample, DataEra);
2121
}
22-
23-
// Setup weight function map for SystematicHelper
24-
//if (!IsDATA && run_syst) {
25-
// setupWeightFunctions();
26-
// }
2722
}
2823

2924
void DiLepton::executeEvent() {
3025
Event ev = GetEvent();
26+
27+
// Initial cutflow entry
28+
float initialWeight = IsDATA ? 1.0 : MCweight() * ev.GetTriggerLumi("Full");
29+
fillCutflow(CutStage::Initial, Channel::NONE, initialWeight, "Central");
30+
3131
RVec<Jet> rawJets = GetAllJets();
3232
if (!PassNoiseFilter(rawJets, ev)) return;
33+
fillCutflow(CutStage::NoiseFilter, Channel::NONE, initialWeight, "Central");
3334

3435
RVec<Muon> rawMuons = GetAllMuons();
35-
if (!RunNoVetoMap && !PassVetoMap(rawJets, rawMuons, "jetvetomap")) return; // For Run3, reject events if any jet is within the veto map
36+
if (!(RunNoVetoMap || PassVetoMap(rawJets, rawMuons, "jetvetomap"))) return; // For Run3, reject events if any jet is within the veto map
37+
fillCutflow(CutStage::VetoMap, Channel::NONE, initialWeight, "Central");
3638

3739
RVec<Electron> rawElectrons = GetAllElectrons();
3840
Particle METv = ev.GetMETVector(Event::MET_Type::PUPPI);
@@ -43,9 +45,12 @@ void DiLepton::executeEvent() {
4345
if (!IsDATA && run_syst && systHelper) {
4446
// Step 1: Process Central objects and weight-only systematics
4547
RecoObjects centralObjects = defineObjects(ev, rawMuons, rawElectrons, rawJets, genJets, METv, "Central");
46-
Channel selectedChannel = selectEvent(ev, centralObjects);
48+
fillCutflow(CutStage::LeptonSelection, Channel::NONE, initialWeight, "Central");
49+
50+
Channel selectedChannel = selectEvent(ev, centralObjects, "Central");
4751

4852
if (selectedChannel != Channel::NONE) {
53+
fillCutflow(CutStage::Final, selectedChannel, initialWeight, "Central");
4954
// Fill Central with nominal weights
5055
WeightInfo centralWeights = getWeights(selectedChannel, ev, centralObjects, genParts, "Central");
5156
fillObjects(selectedChannel, centralObjects, centralWeights, "Central");
@@ -62,7 +67,7 @@ void DiLepton::executeEvent() {
6267

6368
// Define objects with systematic variation
6469
RecoObjects recoObjects = defineObjects(ev, rawMuons, rawElectrons, rawJets, genJets, METv, systName);
65-
Channel systChannel = selectEvent(ev, recoObjects);
70+
Channel systChannel = selectEvent(ev, recoObjects, systName);
6671

6772
if (systChannel != Channel::NONE) {
6873
WeightInfo weights = getWeights(systChannel, ev, recoObjects, genParts, systName);
@@ -73,9 +78,12 @@ void DiLepton::executeEvent() {
7378
} else {
7479
// Process only Central for DATA or when systematics are off
7580
RecoObjects recoObjects = defineObjects(ev, rawMuons, rawElectrons, rawJets, genJets, METv, "Central");
76-
Channel selectedChannel = selectEvent(ev, recoObjects);
81+
fillCutflow(CutStage::LeptonSelection, Channel::NONE, initialWeight, "Central");
82+
83+
Channel selectedChannel = selectEvent(ev, recoObjects, "Central");
7784

7885
if (selectedChannel != Channel::NONE) {
86+
fillCutflow(CutStage::Final, selectedChannel, initialWeight, "Central");
7987
WeightInfo weights = getWeights(selectedChannel, ev, recoObjects, genParts, "Central");
8088
fillObjects(selectedChannel, recoObjects, weights, "Central");
8189
}
@@ -162,48 +170,62 @@ DiLepton::RecoObjects DiLepton::defineObjects(Event& ev,
162170
return objects;
163171
}
164172

165-
DiLepton::Channel DiLepton::selectEvent(Event& ev, const RecoObjects& recoObjects) {
173+
DiLepton::Channel DiLepton::selectEvent(Event& ev, const RecoObjects& recoObjects, const TString& syst) {
166174
const RVec<Muon>& looseMuons = recoObjects.looseMuons;
167175
const RVec<Muon>& tightMuons = recoObjects.tightMuons;
168176
const RVec<Electron>& looseElectrons = recoObjects.looseElectrons;
169177
const RVec<Electron>& tightElectrons = recoObjects.tightElectrons;
170178
const RVec<Jet>& jets = recoObjects.tightJets_vetoLep;
171179
const RVec<Jet>& bjets = recoObjects.bjets;
172180

181+
float weight = IsDATA ? 1.0 : MCweight() * ev.GetTriggerLumi("Full");
173182
bool isDiMu = (tightMuons.size() == 2 && looseMuons.size() == 2 &&
174183
tightElectrons.size() == 0 && looseElectrons.size() == 0);
175184
bool isEMu = (tightMuons.size() == 1 && looseMuons.size() == 1 &&
176185
tightElectrons.size() == 1 && looseElectrons.size() == 1);
177186

178187
if (channel == Channel::DIMU) {
179188
if (!isDiMu) return Channel::NONE;
189+
fillCutflow(CutStage::LeptonSelection, Channel::DIMU, weight, syst);
180190
}
181191
if (channel == Channel::EMU) {
182192
if (!isEMu) return Channel::NONE;
193+
fillCutflow(CutStage::LeptonSelection, Channel::EMU, weight, syst);
183194
}
184195

185196
// DiMu selection
186197
if (channel == Channel::DIMU) {
187198
if (!ev.PassTrigger(DblMuTriggers)) return Channel::NONE;
199+
fillCutflow(CutStage::Trigger, Channel::DIMU, weight, syst);
200+
188201
const Muon& mu1 = tightMuons[0];
189202
const Muon& mu2 = tightMuons[1];
190203
if (mu1.Pt() <= 20.) return Channel::NONE;
191204
if (mu2.Pt() <= 10.) return Channel::NONE;
192205
if (mu1.Charge() + mu2.Charge() != 0) return Channel::NONE;
193206
Particle pair = mu1 + mu2;
194207
if (pair.M() <= 50.) return Channel::NONE;
208+
fillCutflow(CutStage::KinematicCuts, Channel::DIMU, weight, syst);
195209
return Channel::DIMU;
196210
}
197211
// EMu selection
198212
else if (channel == Channel::EMU) {
199213
if (!ev.PassTrigger(EMuTriggers)) return Channel::NONE;
214+
fillCutflow(CutStage::Trigger, Channel::EMU, weight, syst);
215+
200216
const Muon& mu = tightMuons[0];
201217
const Electron& ele = tightElectrons[0];
202218
if (!((mu.Pt() > 20. && ele.Pt() > 15.) || (mu.Pt() > 10. && ele.Pt() > 25.))) return Channel::NONE;
203219
if (mu.Charge() + ele.Charge() != 0) return Channel::NONE;
204220
if (mu.DeltaR(ele) <= 0.4) return Channel::NONE;
221+
fillCutflow(CutStage::KinematicCuts, Channel::EMU, weight, syst);
222+
205223
if (jets.size() < 2) return Channel::NONE;
224+
fillCutflow(CutStage::JetRequirements, Channel::EMU, weight, syst);
225+
206226
if (bjets.size() < 1) return Channel::NONE;
227+
fillCutflow(CutStage::BjetRequirements, Channel::EMU, weight, syst);
228+
207229
return Channel::EMU;
208230
}
209231

@@ -458,3 +480,12 @@ void DiLepton::fillObjects(const DiLepton::Channel& channel, const RecoObjects&
458480
FillHist(Form("%s/%s/pair/mass", channelStr.Data(), syst.Data()), pair.M(), weight, 300, 0., 300.);
459481
}
460482
}
483+
484+
void DiLepton::fillCutflow(CutStage stage, const Channel& channel, float weight, const TString& syst) {
485+
if (syst != "Central") return;
486+
TString channelStr = channelToString(channel);
487+
if (channelStr == "NONE") channelStr = "ALL";
488+
489+
int cutIndex = static_cast<int>(stage);
490+
FillHist(Form("%s/%s/cutflow", channelStr.Data(), syst.Data()), cutIndex, weight, 9, 0., 9.);
491+
}

‎Analyzers/src/DiLeptonBase.cc‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -17,11 +17,11 @@ void DiLeptonBase::initializeAnalyzer() {
1717

1818
// Lepton IDs and triggers
1919
MuonIDs = new IDContainer("HcToWATight", "HcToWALoose");
20-
ElectronIDs = new IDContainer("HcToWATight", "HcToWALoose");
20+
ElectronIDs = new IDContainer("HcToWATight", ((Run == 2) ? "HcToWALooseRun2" : "HcToWALooseRun3"));
2121
if (DataEra == "2016preVFP") {
2222
DblMuTriggers = {
2323
"HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL",
24-
"HLT_Mu17_TrkIsoVVL_TkMu8_TrkIsoVVL "
24+
"HLT_Mu17_TrkIsoVVL_TkMu8_TrkIsoVVL"
2525
};
2626
EMuTriggers = {
2727
"HLT_Mu8_TrkIsoVVL_Ele23_CaloIdL_TrackIdL_IsoVL",

‎Analyzers/src/SKNanoLoader.cc‎

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -895,12 +895,6 @@ void SKNanoLoader::Init() {
895895
SafeSetBranchAddress("Electron_mvaNoIso", Electron_mvaNoIso.data());
896896
SafeSetBranchAddress("Electron_mvaNoIso_WP80", Electron_mvaNoIso_WP80.data());
897897
SafeSetBranchAddress("Electron_mvaNoIso_WP90", Electron_mvaNoIso_WP90.data());
898-
} else if(Run == 2) {
899-
SafeSetBranchAddress("Electron_mvaIso_WPL", Electron_mvaIso_WPL.data());
900-
SafeSetBranchAddress("Electron_mvaNoIso", Electron_mvaNoIso.data());
901-
SafeSetBranchAddress("Electron_mvaNoIso_WP80", Electron_mvaNoIso_WP80.data());
902-
SafeSetBranchAddress("Electron_mvaNoIso_WP90", Electron_mvaNoIso_WP90.data());
903-
SafeSetBranchAddress("Electron_mvaNoIso_WPLoose", Electron_mvaNoIso_WPL.data());
904898
} else if(Run == 2) {
905899
SafeSetBranchAddress("Electron_cutBased", Electron_cutBased_RunII.data());
906900
SafeSetBranchAddress("Electron_genPartIdx", Electron_genPartIdx_RunII.data());

‎DataFormats/include/Electron.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -139,7 +139,8 @@ class Electron : public Lepton {
139139
// Private IDs
140140
bool Pass_CaloIdL_TrackIdL_IsoVL() const;
141141
bool Pass_HcToWABaseline() const;
142-
bool Pass_HcToWALoose() const;
142+
bool Pass_HcToWALooseRun2() const;
143+
bool Pass_HcToWALooseRun3() const;
143144
bool Pass_HcToWATight() const;
144145
// MVA scores
145146
enum class MVATYPE {NONE, MVAISO, MVANOISO, MVATTH};

‎DataFormats/src/Electron.cc‎

Lines changed: 41 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -76,20 +76,21 @@ bool Electron::PassID(const TString ID) const {
7676
if (etaRegion() == ETAREGION::GAP) return false;
7777

7878
// POG
79-
if (ID == "") return true;
80-
if (ID == "NOCUT") return true;
81-
if (ID == "POGVeto") return (int)(CutBased()) >= (int)(WORKINGPOINT::VETO);
82-
if (ID == "POGLoose") return (int)(CutBased()) >= (int)(WORKINGPOINT::LOOSE);
83-
if (ID == "POGMedium") return (int)(CutBased()) >= (int)(WORKINGPOINT::MEDIUM);
84-
if (ID == "POGTight") return (int)(CutBased()) >= (int)(WORKINGPOINT::TIGHT);
85-
if (ID == "POGHEEP") return isCutBasedHEEP();
86-
if (ID == "POGMVAIsoWP80") return isMVAIsoWP80();
87-
if (ID == "POGMVAIsoWP90") return isMVAIsoWP90();
88-
if (ID == "POGMVANoIsoWP80") return isMVANoIsoWP80();
89-
if (ID == "POGMVANoIsoWP90") return isMVANoIsoWP90();
90-
if (ID == "POGMVANoIsoWPLoose") return isMVANoIsoWPLoose();
91-
if (ID == "HcToWATight") return Pass_HcToWATight();
92-
if (ID == "HcToWALoose") return Pass_HcToWALoose();
79+
if (ID == "") return true;
80+
if (ID == "NOCUT") return true;
81+
if (ID == "POGVeto") return (int)(CutBased()) >= (int)(WORKINGPOINT::VETO);
82+
if (ID == "POGLoose") return (int)(CutBased()) >= (int)(WORKINGPOINT::LOOSE);
83+
if (ID == "POGMedium") return (int)(CutBased()) >= (int)(WORKINGPOINT::MEDIUM);
84+
if (ID == "POGTight") return (int)(CutBased()) >= (int)(WORKINGPOINT::TIGHT);
85+
if (ID == "POGHEEP") return isCutBasedHEEP();
86+
if (ID == "POGMVAIsoWP80") return isMVAIsoWP80();
87+
if (ID == "POGMVAIsoWP90") return isMVAIsoWP90();
88+
if (ID == "POGMVANoIsoWP80") return isMVANoIsoWP80();
89+
if (ID == "POGMVANoIsoWP90") return isMVANoIsoWP90();
90+
//if (ID == "POGMVANoIsoWPLoose") return isMVANoIsoWPLoose();
91+
if (ID == "HcToWATight") return Pass_HcToWATight();
92+
if (ID == "HcToWALooseRun2") return Pass_HcToWALooseRun2();
93+
if (ID == "HcToWALooseRun3") return Pass_HcToWALooseRun3();
9394

9495
cerr << "[Electron::PassID] " << ID << " is not implemented" << endl;
9596
exit(ENODATA);
@@ -158,7 +159,7 @@ bool Electron::Pass_HcToWATight() const {
158159
return true;
159160
}
160161

161-
bool Electron::Pass_HcToWALoose() const {
162+
bool Electron::Pass_HcToWALooseRun2() const {
162163
if (! Pass_HcToWABaseline()) return false;
163164
if (! (SIP3D() < 8.)) return false;
164165
if (! (MiniPFRelIso() < 0.4)) return false;
@@ -178,4 +179,28 @@ bool Electron::Pass_HcToWALoose() const {
178179
}
179180
if (! (isMVANoIsoWP90() || passMVAIDNoIsoCut)) return false;
180181
return true;
181-
}
182+
}
183+
184+
bool Electron::Pass_HcToWALooseRun3() const {
185+
if (! Pass_HcToWABaseline()) return false;
186+
if (! (SIP3D() < 8.)) return false;
187+
if (! (MiniPFRelIso() < 0.4)) return false;
188+
const float cutIB=0.5, cutOB=-0.8, cutEC=-0.5;
189+
bool passMVAIDNoIsoCut = false;
190+
switch(etaRegion()) {
191+
case ETAREGION::IB:
192+
if (! (MvaNoIso() > cutIB)) passMVAIDNoIsoCut = true;
193+
break;
194+
case ETAREGION::OB:
195+
if (! (MvaNoIso() > cutOB)) passMVAIDNoIsoCut = true;
196+
break;
197+
case ETAREGION::EC:
198+
if (! (MvaNoIso() > cutEC)) passMVAIDNoIsoCut = true;
199+
break;
200+
default: break;
201+
}
202+
if (! (isMVANoIsoWP90() || passMVAIDNoIsoCut)) return false;
203+
return true;
204+
}
205+
206+

0 commit comments

Comments
 (0)