diff --git a/include/AnalysisOutput.h b/include/AnalysisOutput.h new file mode 100755 index 0000000..962257f --- /dev/null +++ b/include/AnalysisOutput.h @@ -0,0 +1,19 @@ +#pragma once +#include "UHH2/core/include/Event.h" +#include "UHH2/core/include/AnalysisModule.h" + +#include + +using namespace std; + +class WriteOutput: public uhh2::AnalysisModule{ +public: + + explicit WriteOutput(uhh2::Context&); + virtual bool process(uhh2::Event & ) override; + +private: + uhh2::Event::Handleh_weight; + + +}; diff --git a/include/GenSelections.h b/include/GenSelections.h index 38f7fc9..3f72db5 100644 --- a/include/GenSelections.h +++ b/include/GenSelections.h @@ -1,21 +1,9 @@ #pragma once #include -#include #include -#include -#include -#include -#include -#include #include -#include - -#include -#include -#include -#include using namespace std; diff --git a/include/MTopJetUtils.h b/include/MTopJetUtils.h new file mode 100644 index 0000000..83af488 --- /dev/null +++ b/include/MTopJetUtils.h @@ -0,0 +1,9 @@ +#pragma once + +#include + +//// + +const Particle* leading_lepton(const uhh2::Event&); + +//// \ No newline at end of file diff --git a/include/RecoSelections.h b/include/RecoSelections.h index 2cf052f..91b7e20 100644 --- a/include/RecoSelections.h +++ b/include/RecoSelections.h @@ -12,6 +12,9 @@ #include #include #include +#include + +#include #include #include @@ -35,4 +38,51 @@ namespace uhh2 { private: float min_met_, max_met_; }; + + //////////////////////////////////////////////////////////////// + + class ElectronEtaVeto : public Selection { + + public: + explicit ElectronEtaVeto(double, double); + virtual bool passes(const Event&) override; + + private: + double lower, upper; + }; + + //////////////////////////////////////////////////////////////// + + class TwoDCut : public Selection { + + public: + explicit TwoDCut(float min_deltaR, float min_pTrel): min_deltaR_(min_deltaR), min_pTrel_(min_pTrel) {} + virtual bool passes(const Event&) override; + + private: + float min_deltaR_, min_pTrel_; + }; + + //////////////////////////////////////////////////////////////// + + class BadHCALSelection: public uhh2::Selection { + public: + BadHCALSelection(uhh2::Context &ctx, long int seed = 123456789); + virtual bool passes(const uhh2::Event &event) override; + + private: + TRandom *m_rng; + long int m_seed; + Year year; + int m_runnumber = 319077; + double m_lumi_ratio = 0.64844705699; // (Run 319077(17.370008 pb-1) + Run C + Run D) / all 2018 + + double m_interval_eta = -1.3; + double m_interval_phi_low = -1.57; + double m_interval_phi_high = -0.87; + + }; + + //////////////////////////////////////////////////////////////// + } diff --git a/include/RemoveLepton.h b/include/RemoveLepton.h new file mode 100755 index 0000000..1e86b7b --- /dev/null +++ b/include/RemoveLepton.h @@ -0,0 +1,31 @@ +#pragma once +#include "UHH2/core/include/Event.h" +#include "UHH2/common/include/TTbarGen.h" +#include "UHH2/core/include/GenTopJet.h" + +#include + +using namespace std; +using namespace uhh2; + +class RemoveLepton: public uhh2::AnalysisModule{ + + public: + RemoveLepton(uhh2::Context & ctx, const std::string & name_jet); + virtual bool process(uhh2::Event & ) override; + + private: + uhh2::Event::Handle>h_topjets; +}; + + +class RemoveLeptonGen: public uhh2::AnalysisModule{ + + public: + RemoveLeptonGen(uhh2::Context & ctx, const std::string & name_jet); + virtual bool process(uhh2::Event & ) override; + + private: + uhh2::Event::Handle>h_topjets; + uhh2::Event::Handle h_ttbargen; +}; diff --git a/src/AnalysisOutput.cxx b/src/AnalysisOutput.cxx new file mode 100755 index 0000000..6836673 --- /dev/null +++ b/src/AnalysisOutput.cxx @@ -0,0 +1,15 @@ +#include + +using namespace std; + +WriteOutput::WriteOutput(uhh2::Context & ctx): + h_weight(ctx.declare_event_output("weight")) {} + +bool WriteOutput::process(uhh2::Event & event){ + + double weight = event.weight; + + event.set(h_weight, weight); + + return true; +} diff --git a/src/GenSelections.cxx b/src/GenSelections.cxx index 1c3423c..a1d149a 100644 --- a/src/GenSelections.cxx +++ b/src/GenSelections.cxx @@ -1,5 +1,4 @@ #include -#include diff --git a/src/MTopJetPreSelection.cxx b/src/MTopJetPreSelection.cxx index 2719168..54bb274 100644 --- a/src/MTopJetPreSelection.cxx +++ b/src/MTopJetPreSelection.cxx @@ -4,6 +4,7 @@ #include #include #include +#include #include #include @@ -15,19 +16,19 @@ #include #include #include -#include "UHH2/common/include/YearRunSwitchers.h" - -// Hists -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include + // #include #include #include +#include +#include using namespace std; @@ -38,17 +39,39 @@ class MTopJetPreSelection : public ModuleBASE { virtual bool process(uhh2::Event&) override; protected: + enum lepton { muon, elec }; + lepton channel_; + + // cleaners & Correctors + std::unique_ptr muoSR_cleaner; + std::unique_ptr eleSR_cleaner; + std::unique_ptr common; + std::unique_ptr jet_cleaner1; + std::unique_ptr jet_cleaner2; // selections - std::unique_ptr lumi_sel; + std::unique_ptr remove_lepton_rec; + std::unique_ptr remove_lepton_gen; + std::unique_ptr lumi_sel; std::unique_ptr genmttbar_sel; std::unique_ptr met_sel; std::unique_ptr muon_sel; std::unique_ptr elec_sel; + std::unique_ptr elec_etaveto; + std::unique_ptr elec_sel_120; + std::unique_ptr elec_sel_triggerA; std::unique_ptr SemiLepDecay; std::unique_ptr GenMuonPT; std::unique_ptr GenElecPT; + std::unique_ptr pv_sel; + std::unique_ptr twodcut_sel; + std::unique_ptr badhcal_sel; + std::unique_ptr trigger_mu_A; + std::unique_ptr trigger_mu_B; + std::unique_ptr trigger_el_A; + std::unique_ptr trigger_el_B; + std::unique_ptr trigger_el_C; std::unique_ptr ttgenprod; @@ -56,10 +79,24 @@ class MTopJetPreSelection : public ModuleBASE { Event::Handleh_recsel; Event::Handleh_gensel; Event::Handle>h_fatjets; + Event::Handle>h_genfatjets; + + //write output + std::unique_ptr output; + + Year year; + uint nJets; // bools bool isMC; + bool isTTbar; + bool isElectronStream; + bool isPhotonStream; bool debug; + + bool year_16; + bool year_17; + bool year_18; }; @@ -74,12 +111,50 @@ class MTopJetPreSelection : public ModuleBASE { MTopJetPreSelection::MTopJetPreSelection(uhh2::Context& ctx){ debug = string2bool(ctx.get("Debug","false")); // look for Debug, expect false if not found + nJets = atoi(ctx.get("nJets","2").c_str()); // look for nJets, expect 2 if not found + + //// YEARSWITCHER + year = extract_year(ctx); // Ask for the year of Event + + cout << "post extract_year" << endl; + + year_16 = false; + year_17 = false; + year_18 = false; + if(year == Year::is2016v3) year_16 = true; + else if(year == Year::isUL17) year_17 = true; + else if(year == Year::isUL18) year_18 = true; + else throw runtime_error("In MTopJetPreSelection: This Event is not from 2016v3, 2017v2 or 2018!"); //// CONFIGURATION if(debug) cout << "Configuration" << endl; isMC = (ctx.get("dataset_type") == "MC"); + TString dataset_version = (TString) ctx.get("dataset_version"); + if(dataset_version.Contains("TTbar") || dataset_version.Contains("TTTo")) isTTbar = true; + else isTTbar = false; + + if(dataset_version.Contains("SingleElec")) isElectronStream = true; + else isElectronStream = false; + + if(dataset_version.Contains("SinglePhoton")) isPhotonStream = true; + else isPhotonStream = false; + + if(isElectronStream) cout << "ELECTRON STREAM DETECTED" << endl; + if(isPhotonStream) cout << "PHOTON STREAM DETECTED" << endl; + + const std::string& channel = ctx.get("channel", ""); + if (channel == "muon") channel_ = muon; + else if(channel == "elec") channel_ = elec; + else { + + std::string log("TTbarLJAnalysisLiteModule::TTbarLJAnalysisLiteModule -- "); + log += "invalid argument for 'channel' key in xml file (must be 'muon' or 'elec'): \""+channel+"\""; + + throw std::runtime_error(log); + } + //// HANDLES if(debug) cout << "Output and Handles" << endl; @@ -87,9 +162,13 @@ MTopJetPreSelection::MTopJetPreSelection(uhh2::Context& ctx){ h_gensel = ctx.declare_event_output("passed_gensel"); h_fatjets = ctx.get_handle>("xconeCHS"); + h_genfatjets=ctx.get_handle>("genXCone33TopJets"); + if(nJets == 3) { + //unfortunate naming, but genXCone3TopJets has 3 clusterd fatjets, while genXCone33TopJets has 2 fatjets + h_genfatjets=ctx.get_handle>("genXCone3TopJets"); + } //// COMMON MODULES - if(debug) cout << "Common Modules" << endl; if(!isMC) lumi_sel.reset(new LumiSelection(ctx)); @@ -109,10 +188,59 @@ MTopJetPreSelection::MTopJetPreSelection(uhh2::Context& ctx){ genmttbar_sel.reset(new uhh2::AndSelection(ctx)); } + //// IDs + MuonId muid = AndId(MuonID(Muon::Tight), PtEtaCut(55., 2.4)); + ElectronId eleid_noiso55; + if(year_16) eleid_noiso55 = AndId(PtEtaSCCut(55., 2.4), ElectronTagID(Electron::mvaEleID_Spring16_GeneralPurpose_V1_wp90)); + else eleid_noiso55 = AndId(PtEtaSCCut(55., 2.4), ElectronTagID(Electron::mvaEleID_Fall17_noIso_V2_wp90)); + ElectronId eleid_noiso120; + if(year_16) eleid_noiso120 = AndId(PtEtaSCCut(120., 2.4), ElectronTagID(Electron::mvaEleID_Spring16_GeneralPurpose_V1_wp90)); + else eleid_noiso120 = AndId(PtEtaSCCut(120., 2.4), ElectronTagID(Electron::mvaEleID_Fall17_noIso_V2_wp90)); + ElectronId eleid_iso55; + if(year_16) eleid_iso55 = AndId(PtEtaSCCut(55., 2.4), ElectronTagID(Electron::mvaEleID_Spring16_GeneralPurpose_V1_wp90)); + else eleid_iso55 = AndId(PtEtaSCCut(55., 2.4), ElectronTagID(Electron::mvaEleID_Fall17_iso_V2_wp90)); + JetId jetid_cleaner = AndId(JetPFID(JetPFID::WP_TIGHT_CHS), PtEtaCut(30.0, 2.4)); + + //// CLEANER + muoSR_cleaner.reset(new MuonCleaner(muid)); + eleSR_cleaner.reset(new ElectronCleaner(eleid_noiso55)); + + common.reset(new CommonModules()); + common->set_HTjetid(jetid_cleaner); + common->switch_jetlepcleaner(true); + common->switch_metcorrection(); + common->disable_mcpileupreweight(); + common->init(ctx); + + jet_cleaner1.reset(new JetCleaner(ctx, 15., 3.0)); + jet_cleaner2.reset(new JetCleaner(ctx, 30., 2.4)); + //// EVENT SELECTION REC - met_sel.reset(new METCut(40, uhh2::infinity)); - muon_sel.reset(new NMuonSelection(1, -1, MuonId(PtEtaCut(50, 2.4 )))); - elec_sel.reset(new NElectronSelection(1, -1, ElectronId(PtEtaCut(50, 2.4)))); + met_sel.reset(new METCut(50 , uhh2::infinity)); + if(channel_ == elec){ + muon_sel.reset(new NMuonSelection(0, 0, muid)); + elec_sel.reset(new NElectronSelection(1, 1, eleid_noiso55)); + } + else if (channel_ == muon){ + muon_sel.reset(new NMuonSelection(1, 1, muid)); + elec_sel.reset(new NElectronSelection(0, 0, eleid_noiso55)); + } + elec_etaveto.reset(new ElectronEtaVeto(1.44, 1.57)); + pv_sel.reset(new NPVSelection(1, -1, PrimaryVertexId(StandardPrimaryVertexId()))); + elec_sel_triggerA.reset(new NElectronSelection(1, 1, eleid_iso55)); + elec_sel_120.reset(new NElectronSelection(1, 1, eleid_noiso120)); + twodcut_sel.reset(new TwoDCut(0.4, 40)); + badhcal_sel.reset(new BadHCALSelection(ctx)); + + //// TRIGGER + trigger_mu_A = uhh2::make_unique("HLT_Mu50_v*"); + trigger_mu_B = uhh2::make_unique("HLT_TkMu50_v*"); + if(year_16) trigger_el_A = uhh2::make_unique("HLT_Ele27_WPTight_Gsf_v*"); + else if(year_17) trigger_el_A = uhh2::make_unique("HLT_Ele35_WPTight_Gsf_v*"); + else if(year_18) trigger_el_A = uhh2::make_unique("HLT_Ele32_WPTight_Gsf_v*"); + trigger_el_B = uhh2::make_unique("HLT_Ele115_CaloIdVT_GsfTrkIdT_v*"); + if(year_16) trigger_el_C = uhh2::make_unique("HLT_Photon175_v*"); + else trigger_el_C = uhh2::make_unique("HLT_Photon200_v*"); //// EVENTS SELECTION GEN SemiLepDecay.reset(new TTbarSemilep(ctx)); @@ -121,6 +249,14 @@ MTopJetPreSelection::MTopJetPreSelection(uhh2::Context& ctx){ GenElecPT.reset(new GenElecSel(ctx, 55.)); } + //// PRODUCER + remove_lepton_rec.reset(new RemoveLepton(ctx, "xconeCHS")); + if(isMC) { + if(nJets == 2) remove_lepton_gen.reset(new RemoveLeptonGen(ctx, "genXCone33TopJets")); + else if(nJets == 3) remove_lepton_gen.reset(new RemoveLeptonGen(ctx, "genXCone3TopJets")); + } + output.reset(new WriteOutput(ctx)); + } /* @@ -135,33 +271,171 @@ bool MTopJetPreSelection::process(uhh2::Event& event){ if(debug) cout << "Start Module - Process" << endl; + // get Rec and Gen jets + std::vector jets = event.get(h_fatjets); + std::vector genJets = event.get(h_genfatjets); + + // check number of jets + if(jets.size() < nJets) return false; + if(!event.isRealData){ + if(genJets.size() < nJets) return false; + } + + if(debug) cout << "luminosity sections and GEN M-ttbar selection" << endl; + /* CMS-certified luminosity sections */ if(event.isRealData){ if(!lumi_sel->passes(event)) return false; } + /* GEN M-ttbar selection */ if(!event.isRealData){ - /* GEN M-ttbar selection */ ttgenprod->process(event); if(!genmttbar_sel->passes(event)) return false; } + //// CLEANER + muoSR_cleaner->process(event); + if(event.muons->size() > 0) sort_by_pt(*event.muons); + + eleSR_cleaner->process(event); + if(event.electrons->size() > 0) sort_by_pt(*event.electrons); + + if(!common->process(event)) return false; + + jet_cleaner1->process(event); + sort_by_pt(*event.jets); + //// EVENT SELECTION REC bool passed_recsel; bool pass_lep_number = ((event.muons->size() >= 1) || (event.electrons->size() >= 1)); - bool pass_lepsel = (muon_sel->passes(event) || elec_sel->passes(event)); + bool pass_lepsel = (muon_sel->passes(event) && elec_sel->passes(event)); + if(channel_ == elec) pass_lepsel = elec_etaveto->passes(event); bool pass_met = met_sel->passes(event); + bool pass_prim_vert = pv_sel->passes(event); // cut on fatjet pt - bool passed_fatpt=false; - std::vector jets = event.get(h_fatjets); + bool pass_fatpt=false; double ptcut = 200; for(auto jet: jets){ - if(jet.pt() > ptcut) passed_fatpt = true; + if(jet.pt() > ptcut) pass_fatpt = true; } - if(pass_lep_number && pass_met && pass_lepsel && passed_fatpt) passed_recsel = true; + // trigger + bool pass_lep_trigger = false; + bool pass_mu_trigger = false; + bool pass_elec_trigger = true; + + bool elec_is_isolated = false; + + // for DATA until run 274954 -> use only Trigger A + // for MC and DATA from 274954 -> use "A || B" + // HLT_TkMu50_v is not availible for 2017&18. Alternative still needs to be inlcuded !!!!!!!!!!!!!!!! + if(channel_ == muon){ + if(year_16){ + if(!isMC && event.run < 274954){ + pass_mu_trigger = trigger_mu_A->passes(event); + } + else{ + pass_mu_trigger = (trigger_mu_A->passes(event) || trigger_mu_B->passes(event)); + } + } + if(year_17 || year_18){ + pass_mu_trigger = trigger_mu_A->passes(event); + } + } + // only use triggerA and isolation if elec pt < 120 + // for pt > 120 use triggerB || triggerC + else if(channel_ == elec){ + // pt < 120 + if(!elec_sel_120->passes(event)){ + if(isPhotonStream) pass_elec_trigger = false; + if(!trigger_el_A->passes(event)) pass_elec_trigger = false; + if(!elec_sel_triggerA->passes(event)) pass_elec_trigger = false; + if(pass_elec_trigger) elec_is_isolated = true; + } + // pt > 120 + else{ + // MC: A || B + if(isMC){ + if(!trigger_el_B->passes(event) && !trigger_el_C->passes(event)) pass_elec_trigger = false; + } + // DATA 2016: if elec B, if photon !B && C + else{ + if(year_16){ + if(isPhotonStream){ + if(trigger_el_B->passes(event)) pass_elec_trigger = false; + if(!trigger_el_C->passes(event)) pass_elec_trigger = false; + } + else if(isElectronStream){ + if(!trigger_el_B->passes(event)) pass_elec_trigger = false; + } + } + // DATA 2017: if elec B, if photon !B && C + else if(year_17){ + // Trigger B does not exist in 2017B, so use A here + if(event.run <= 299329){ + if(isPhotonStream){ + if(trigger_el_A->passes(event)) pass_elec_trigger = false; + if(!trigger_el_C->passes(event)) pass_elec_trigger = false; + } + else if(isElectronStream){ + if(!trigger_el_A->passes(event)) pass_elec_trigger = false; + if(!elec_sel_triggerA->passes(event)) pass_elec_trigger = false; + if(pass_elec_trigger) elec_is_isolated = true; + } + } + // For all other runs, same as in 2016 + else{ + if(isPhotonStream){ + if(trigger_el_B->passes(event)) pass_elec_trigger = false; + if(!trigger_el_C->passes(event)) pass_elec_trigger = false; + } + else if(isElectronStream){ + if(!trigger_el_B->passes(event)) pass_elec_trigger = false; + } + } + } + // DATA 2018: B || C + else if(year_18){ + if(!trigger_el_B->passes(event) && !trigger_el_C->passes(event)) pass_elec_trigger = false; + } + } + } + } + + pass_lep_trigger = pass_mu_trigger && pass_elec_trigger; + + // lepton-2Dcut variables + bool pass_twodcut = twodcut_sel->passes(event); + + for(auto& muo : *event.muons){ + float dRmin, pTrel; + std::tie(dRmin, pTrel) = drmin_pTrel(muo, *event.jets); + + muo.set_tag(Muon::twodcut_dRmin, dRmin); + muo.set_tag(Muon::twodcut_pTrel, pTrel); + } + + for(auto& ele : *event.electrons){ + + float dRmin, pTrel; + std::tie(dRmin, pTrel) = drmin_pTrel(ele, *event.jets); + + ele.set_tag(Electron::twodcut_dRmin, dRmin); + ele.set_tag(Electron::twodcut_pTrel, pTrel); + } + if(elec_is_isolated) pass_twodcut = true; // do not do 2D cut for isolated electrons + + // second cleaner + jet_cleaner2->process(event); + sort_by_pt(*event.jets); + + // bad HCAL + bool pass_badhcal = badhcal_sel->passes(event); + + if(pass_lep_number && pass_met && pass_lepsel && pass_fatpt && pass_prim_vert && pass_lep_trigger && pass_twodcut && pass_badhcal) passed_recsel = true; else passed_recsel = false; //// EVNET SELECTION GEN @@ -187,6 +461,14 @@ bool MTopJetPreSelection::process(uhh2::Event& event){ if(!passed_recsel && !passed_gensel) return false; + //// PRODUCER + remove_lepton_rec->process(event); + if(isTTbar){ + remove_lepton_gen->process(event); + } + + output->process(event); + event.set(h_recsel, passed_recsel); event.set(h_gensel, passed_gensel); diff --git a/src/MTopJetUtils.cxx b/src/MTopJetUtils.cxx new file mode 100644 index 0000000..f253f05 --- /dev/null +++ b/src/MTopJetUtils.cxx @@ -0,0 +1,26 @@ +#include +#include +#include +#include +#include + +#include + + +//// + +const Particle* leading_lepton(const uhh2::Event& event){ + + const Particle* lep(0); + + float ptL_max(0.); + if(event.muons) { for(const auto& mu : *event.muons) { if(mu.pt() > ptL_max){ ptL_max = mu.pt(); lep = μ } } } + if(event.electrons){ for(const auto& el : *event.electrons){ if(el.pt() > ptL_max){ ptL_max = el.pt(); lep = ⪙ } } } + + if(!lep) throw std::runtime_error("leading_lepton -- pt-leading lepton not found"); + + return lep; +} + +//// + diff --git a/src/RecoSelections.cxx b/src/RecoSelections.cxx index 0a8f23a..6d4dd65 100644 --- a/src/RecoSelections.cxx +++ b/src/RecoSelections.cxx @@ -10,4 +10,64 @@ bool uhh2::METCut::passes(const uhh2::Event& event){ float MET = event.met->pt(); return (MET > min_met_) && (MET < max_met_); } + +//////////////////////////////////////////////////////// + +uhh2::ElectronEtaVeto::ElectronEtaVeto(double lower_, double upper_): +lower(lower_), +upper(upper_) {} + +bool uhh2::ElectronEtaVeto::passes(const uhh2::Event& event){ + bool pass_veto = false; + if(event.electrons->size() != 0){ + double eta = fabs(event.electrons->at(0).eta()); + if(eta < lower || eta > upper) pass_veto = true; + } + return pass_veto; +} + +//////////////////////////////////////////////////////// + +bool uhh2::TwoDCut::passes(const uhh2::Event& event){ + + if(event.muons->size() != 0 || event.electrons->size() != 0){ + assert(event.muons && event.electrons && event.jets); + + const Particle* lepton = leading_lepton(event); + + float drmin, ptrel; + std::tie(drmin, ptrel) = drmin_pTrel(*lepton, *event.jets); + + return ((drmin > min_deltaR_) || (ptrel > min_pTrel_)); + } + else return true; +} + +//////////////////////////////////////////////////////// + +BadHCALSelection::BadHCALSelection(Context &ctx, long int seed) { + m_seed = seed; + m_rng = new TRandomMixMax(); + m_rng->SetSeed(m_seed); + year = extract_year(ctx); +} + +bool BadHCALSelection::passes(const Event &event) { + + if (year != Year::is2018) return true; + + // check if event should be removed: + // for data: if event is affected by HEM15/16 + // for mc: draw random sample according to lumi ratio of affected data + if ((event.isRealData && event.run >= m_runnumber) || (!event.isRealData && m_rng->Uniform() < m_lumi_ratio)){ + for (const Electron & e : *event.electrons){ + if (e.eta() < m_interval_eta && e.phi() > m_interval_phi_low && e.phi() < m_interval_phi_high) return false; + } + for (const Jet & j : *event.jets){ + if (j.eta() < m_interval_eta && j.phi() > m_interval_phi_low && j.phi() < m_interval_phi_high) return false; + } + } + return true; +} + //////////////////////////////////////////////////////// diff --git a/src/RemoveLepton.cxx b/src/RemoveLepton.cxx new file mode 100755 index 0000000..5e9f77e --- /dev/null +++ b/src/RemoveLepton.cxx @@ -0,0 +1,133 @@ +#include + +RemoveLepton::RemoveLepton(uhh2::Context & ctx, const std::string & name_jet): +h_topjets(ctx.get_handle>(name_jet)){} + +bool RemoveLepton::process(uhh2::Event & event){ + std::vector topjets = event.get(h_topjets); + std::vector leptons; + for(unsigned int i=0; isize(); i++){ + leptons.push_back(event.muons->at(i)); + } + for(unsigned int i=0; isize(); i++){ + leptons.push_back(event.electrons->at(i)); + } + if(topjets.size() == 0 || leptons.size() == 0) return true; + + for(unsigned int i=0; i new_subjets = topjets[i].subjets(); + for(unsigned int j=0; j 5 && deltaR(new_v4_top, topjets[i].v4()) > M_PI/2){ + cout << "Warning: subtracting lepton flipped jet direction" << endl; + topjets[i].set_v4(LorentzVector()); + } + else topjets[i].set_v4(new_v4_top); + + // if there is also a subjet within 0.4, find subjet that is closed to + // the lepton and subtract the lepton from this one subjet only + double dRmin = 0.4; + int index = 100; + for(unsigned int k=0; k 5 && deltaR(new_v4_sub, new_subjets[index].v4()) > M_PI/2){ + cout << "Warning: subtracting lepton flipped jet direction" << endl; + new_subjets[index].set_v4(LorentzVector()); + } + else new_subjets[index].set_v4(new_v4_sub); + } + } + } + topjets[i].set_subjets(new_subjets); + } + event.set(h_topjets, topjets); + + return true; +} + +/* +.██████ ███████ ███ ██ +██ ██ ████ ██ +██ ███ █████ ██ ██ ██ +██ ██ ██ ██ ██ ██ +.██████ ███████ ██ ████ +*/ + + +RemoveLeptonGen::RemoveLeptonGen(uhh2::Context & ctx, const std::string & name_jet): +h_topjets(ctx.get_handle>(name_jet)), +h_ttbargen(ctx.get_handle("ttbargen")){} + +bool RemoveLeptonGen::process(uhh2::Event & event){ + std::vector topjets = event.get(h_topjets); + GenParticle lepton; + std::vector leptons; + const auto & ttbargen = event.get(h_ttbargen); + if(ttbargen.IsSemiLeptonicDecay()) { + lepton = ttbargen.ChargedLepton(); + leptons.push_back(lepton); + } + + if(topjets.size() == 0 || leptons.size() == 0) return true; + + // for GenTopJet there is no method 'set_subjets' + // So, you need to create new GenTopJets without subjets and add the + // corrected subjets in the end + std::vector new_topjets; + for(unsigned int i=0; i new_subjets = topjets[i].subjets(); + GenTopJet new_topjet; + new_topjet.set_v4(topjets[i].v4()); + // loop over leptons and subtract v4 from topjet if dR < 1.2 + if(deltaR(new_topjet, lepton) < 1.2){ + LorentzVector new_v4_top = new_topjet.v4() - lepton.v4(); + if(new_v4_top.pt() > 5 && deltaR(new_v4_top, new_topjet.v4()) > M_PI/2){ + cout << "Warning: subtracting lepton flipped jet direction" << endl; + new_topjet.set_v4(LorentzVector()); + } + else new_topjet.set_v4(new_v4_top); + + // if there is also a subjet within 0.4, find subjet that is closed to + // the lepton and subtract the lepton from this one subjet only + double dRmin = 0.4; + int index = 100; + for(unsigned int k=0; k 5 && deltaR(new_v4_sub, new_subjets[index].v4()) > M_PI/2){ + cout << "Warning: subtracting lepton flipped jet direction" << endl; + new_subjets[index].set_v4(LorentzVector()); + } + else new_subjets[index].set_v4(new_v4_sub); + } + } + for(auto subjet: new_subjets) new_topjet.add_subjet(subjet); + + new_topjet.set_tau1(topjets.at(i).tau1()); + new_topjet.set_tau2(topjets.at(i).tau2()); + new_topjet.set_tau3(topjets.at(i).tau3()); + new_topjet.set_tau4(topjets.at(i).tau4()); + new_topjet.set_tau1_groomed(topjets.at(i).tau1_groomed()); + new_topjet.set_tau2_groomed(topjets.at(i).tau2_groomed()); + new_topjet.set_tau3_groomed(topjets.at(i).tau3_groomed()); + new_topjet.set_tau4_groomed(topjets.at(i).tau4_groomed()); + + new_topjets.push_back(new_topjet); + } + + event.set(h_topjets, new_topjets); + return true; +}