diff --git a/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py b/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py index ee2867025e2d6..09c5c8169508b 100644 --- a/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py +++ b/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py @@ -809,5 +809,16 @@ Plot1D('rawMVAoldDMdR032017v2', 'rawMVAoldDMdR032017v2', 20, -1, 1, 'byIsolationMVArun2017v2DBoldDMdR0p3wLT raw output discriminator (2017v2)'), ) ), + L1PreFiringWeight = cms.PSet( + sels = cms.PSet(), + plots = cms.VPSet( + Plot1D('Nom', 'Nom', 21, 0.8, 1.01, 'L1 prefiring weight nominal'), + Plot1D('Up', 'Up', 21, 0.8, 1.01, 'L1 prefiring weight uncertainy up'), + Plot1D('Dn', 'Dn', 21, 0.8, 1.01, 'L1 prefiring weight uncertainty down'), + Plot1D('ECAL_Nom', 'ECAL_Nom', 21, 0.8, 1.01, 'L1 prefiring weight for ECAL objects nominal'), + Plot1D('Muon_Nom', 'Muon_Nom', 21, 0.8, 1.01, 'L1 prefiring weight for muons nominal'), + ) + + ), ) ) diff --git a/PhysicsTools/NanoAOD/python/nano_eras_cff.py b/PhysicsTools/NanoAOD/python/nano_eras_cff.py index cfe39c6f1cc4e..62b2723a2fe31 100644 --- a/PhysicsTools/NanoAOD/python/nano_eras_cff.py +++ b/PhysicsTools/NanoAOD/python/nano_eras_cff.py @@ -4,9 +4,11 @@ from Configuration.Eras.Modifier_run2_jme_2016_cff import run2_jme_2016 from Configuration.Eras.Modifier_run2_jme_2017_cff import run2_jme_2017 from Configuration.Eras.Modifier_run2_muon_2016_cff import run2_muon_2016 +from Configuration.Eras.Modifier_run2_muon_2018_cff import run2_muon_2018 from Configuration.Eras.Modifier_run2_HLTconditions_2016_cff import run2_HLTconditions_2016 from Configuration.Eras.Modifier_run2_HLTconditions_2017_cff import run2_HLTconditions_2017 +from Configuration.Eras.Modifier_run2_HLTconditions_2018_cff import run2_HLTconditions_2018 from Configuration.Eras.Modifier_run2_miniAOD_80XLegacy_cff import run2_miniAOD_80XLegacy from Configuration.Eras.Modifier_run2_nanoAOD_92X_cff import run2_nanoAOD_92X diff --git a/PhysicsTools/NanoAOD/python/triggerObjects_cff.py b/PhysicsTools/NanoAOD/python/triggerObjects_cff.py index 6c56743440c09..b398cab6a20e1 100644 --- a/PhysicsTools/NanoAOD/python/triggerObjects_cff.py +++ b/PhysicsTools/NanoAOD/python/triggerObjects_cff.py @@ -216,24 +216,37 @@ selections = selections2016 ) -from PhysicsTools.PatUtils.L1ECALPrefiringWeightProducer_cff import prefiringweight -#Next line will be updated once we get UL2016 maps -run2_jme_2016.toModify( prefiringweight, DataEra = cms.string("2016BtoH")) +from PhysicsTools.PatUtils.L1PrefiringWeightProducer_cff import prefiringweight +#Next lines are for UL2016 maps +(run2_muon_2016 & tracker_apv_vfp30_2016).toModify( prefiringweight, DataEraECAL = cms.string("UL2016preVFP"), DataEraMuon = cms.string("2016preVFP")) +(run2_muon_2016 & ~tracker_apv_vfp30_2016).toModify( prefiringweight, DataEraECAL = cms.string("UL2016postVFP"), DataEraMuon = cms.string("2016postVFP")) #Next line is for UL2017 maps -run2_jme_2017.toModify( prefiringweight, DataEra = cms.string("UL2017BtoF")) +run2_jme_2017.toModify( prefiringweight, DataEraECAL = cms.string("UL2017BtoF"), DataEraMuon = cms.string("20172018")) +#Next line is for UL2018 maps +run2_muon_2018.toModify( prefiringweight, DataEraECAL = cms.string("None"), DataEraMuon = cms.string("20172018")) + #For pre-UL 2017 reprocessing, one should use the original maps and no muon jet protection for modifier in run2_nanoAOD_94XMiniAODv1, run2_nanoAOD_94XMiniAODv2: - modifier.toModify( prefiringweight, DataEra = cms.string("2017BtoF")) + modifier.toModify( prefiringweight, DataEraECAL = cms.string("2017BtoF"), DataEraMuon = cms.string("20172018")) modifier.toModify( prefiringweight, JetMaxMuonFraction = cms.double(-1.) ) #For pre-UL 2016 reprocessing, same thing -run2_nanoAOD_94X2016.toModify( prefiringweight, DataEra = cms.string("2016BtoH") ) +run2_nanoAOD_94X2016.toModify( prefiringweight, DataEraECAL = cms.string("2016BtoH"), DataEraMuon = cms.string("2016") ) run2_nanoAOD_94X2016.toModify( prefiringweight, JetMaxMuonFraction = cms.double(-1.) ) l1PreFiringEventWeightTable = cms.EDProducer("GlobalVariablesTableProducer", + name = cms.string("L1PreFiringWeight"), variables = cms.PSet( - L1PreFiringWeight_Nom = ExtVar(cms.InputTag("prefiringweight:nonPrefiringProb"), "double", doc = "L1 pre-firing event correction weight (1-probability)", precision=8), - L1PreFiringWeight_Up = ExtVar(cms.InputTag("prefiringweight:nonPrefiringProbUp"), "double", doc = "L1 pre-firing event correction weight (1-probability), up var.", precision=8), - L1PreFiringWeight_Dn = ExtVar(cms.InputTag("prefiringweight:nonPrefiringProbDown"), "double", doc = "L1 pre-firing event correction weight (1-probability), down var.", precision=8), + Nom = ExtVar(cms.InputTag("prefiringweight:nonPrefiringProb"), "float", doc = "L1 pre-firing event correction weight (1-probability)", precision=8), + Up = ExtVar(cms.InputTag("prefiringweight:nonPrefiringProbUp"), "float", doc = "L1 pre-firing event correction weight (1-probability), up var.", precision=8), + Dn = ExtVar(cms.InputTag("prefiringweight:nonPrefiringProbDown"), "float", doc = "L1 pre-firing event correction weight (1-probability), down var.", precision=8), + Muon_Nom = ExtVar(cms.InputTag("prefiringweight:nonPrefiringProbMuon"), "float", doc = "Muon L1 pre-firing event correction weight (1-probability)", precision=8), + Muon_SystUp = ExtVar(cms.InputTag("prefiringweight:nonPrefiringProbMuonSystUp"), "float", doc = "Muon L1 pre-firing event correction weight (1-probability), up var. syst.", precision=8), + Muon_SystDn = ExtVar(cms.InputTag("prefiringweight:nonPrefiringProbMuonSystDown"), "float", doc = "Muon L1 pre-firing event correction weight (1-probability), down var. syst.", precision=8), + Muon_StatUp = ExtVar(cms.InputTag("prefiringweight:nonPrefiringProbMuonStatUp"), "float", doc = "Muon L1 pre-firing event correction weight (1-probability), up var. stat.", precision=8), + Muon_StatDn = ExtVar(cms.InputTag("prefiringweight:nonPrefiringProbMuonStatDown"), "float", doc = "Muon L1 pre-firing event correction weight (1-probability), down var. stat.", precision=8), + ECAL_Nom = ExtVar(cms.InputTag("prefiringweight:nonPrefiringProbECAL"), "float", doc = "ECAL L1 pre-firing event correction weight (1-probability)", precision=8), + ECAL_Up = ExtVar(cms.InputTag("prefiringweight:nonPrefiringProbECALUp"), "float", doc = "ECAL L1 pre-firing event correction weight (1-probability), up var.", precision=8), + ECAL_Dn = ExtVar(cms.InputTag("prefiringweight:nonPrefiringProbECALDown"), "float", doc = "ECAL L1 pre-firing event correction weight (1-probability), down var.", precision=8), ) ) @@ -241,4 +254,4 @@ _triggerObjectTables_withL1PreFiring = triggerObjectTables.copy() _triggerObjectTables_withL1PreFiring.replace(triggerObjectTable, prefiringweight + l1PreFiringEventWeightTable + triggerObjectTable) -(run2_HLTconditions_2016 | run2_HLTconditions_2017).toReplaceWith(triggerObjectTables, _triggerObjectTables_withL1PreFiring) +(run2_HLTconditions_2016 | run2_HLTconditions_2017 | run2_HLTconditions_2018).toReplaceWith(triggerObjectTables, _triggerObjectTables_withL1PreFiring) diff --git a/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py b/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py index 8b04ba50a66f9..80614e2c9d1e2 100644 --- a/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py +++ b/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py @@ -537,13 +537,17 @@ def _add_deepMET(process): process.load("RecoEgamma.EgammaTools.slimmedEgammaFromMultiCl_cff") phase2_hgcal.toModify(task, func=lambda t: t.add(process.slimmedEgammaFromMultiClTask)) - # L1 pre-firing weights for 2016 and 2017 + # L1 pre-firing weights for 2016, 2017, and 2018 from Configuration.Eras.Modifier_run2_L1prefiring_cff import run2_L1prefiring - from Configuration.Eras.Modifier_stage1L1Trigger_cff import stage1L1Trigger + from Configuration.Eras.Modifier_stage2L1Trigger_cff import stage2L1Trigger from Configuration.Eras.Modifier_stage2L1Trigger_2017_cff import stage2L1Trigger_2017 - process.load("PhysicsTools.PatUtils.L1ECALPrefiringWeightProducer_cff") - stage1L1Trigger.toModify(process.prefiringweight, DataEra = "2016BtoH") - stage2L1Trigger_2017.toModify(process.prefiringweight, DataEra = "2017BtoF") + from Configuration.Eras.Modifier_stage2L1Trigger_2018_cff import stage2L1Trigger_2018 + from Configuration.Eras.Modifier_tracker_apv_vfp30_2016_cff import tracker_apv_vfp30_2016 + process.load("PhysicsTools.PatUtils.L1PrefiringWeightProducer_cff") + (stage2L1Trigger & tracker_apv_vfp30_2016).toModify(process.prefiringweight, DataEraECAL = "UL2016preVFP", DataEraMuon = "2016preVFP" ) + (stage2L1Trigger & ~tracker_apv_vfp30_2016).toModify(process.prefiringweight, DataEraECAL = "UL2016postVFP", DataEraMuon = "2016postVFP" ) + stage2L1Trigger_2017.toModify(process.prefiringweight, DataEraECAL = "UL2017BtoF", DataEraMuon = "20172018") + stage2L1Trigger_2018.toModify(process.prefiringweight, DataEraECAL = "None", DataEraMuon = "20172018") run2_L1prefiring.toModify(task, func=lambda t: t.add(process.prefiringweight)) from PhysicsTools.PatAlgos.producersHeavyIons.heavyIonJetSetup import removeL1FastJetJECs diff --git a/PhysicsTools/PatUtils/plugins/L1ECALPrefiringWeightProducer.cc b/PhysicsTools/PatUtils/plugins/L1ECALPrefiringWeightProducer.cc deleted file mode 100644 index b67c0c4837326..0000000000000 --- a/PhysicsTools/PatUtils/plugins/L1ECALPrefiringWeightProducer.cc +++ /dev/null @@ -1,231 +0,0 @@ -// -*- C++ -*- -// -// Package: ProdTutorial/L1ECALPrefiringWeightProducer -// Class: L1ECALPrefiringWeightProducer -// -/**\class L1ECALPrefiringWeightProducer L1ECALPrefiringWeightProducer.cc ProdTutorial/L1ECALPrefiringWeightProducer/plugins/L1ECALPrefiringWeightProducer.cc - - Description: [one line class summary] - - Implementation: - [Notes on implementation] -*/ -// -// Original Author: localusers user -// Created: Thu, 08 Nov 2018 16:16:00 GMT -// -// - -// system include files -#include - -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/global/EDProducer.h" - -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/MakerMacros.h" - -#include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "FWCore/Utilities/interface/StreamID.h" - -#include "DataFormats/PatCandidates/interface/Jet.h" -#include "DataFormats/PatCandidates/interface/Photon.h" -#include "TH2.h" - -#include -enum fluctuations { central = 0, up, down }; - -class L1ECALPrefiringWeightProducer : public edm::global::EDProducer<> { -public: - explicit L1ECALPrefiringWeightProducer(const edm::ParameterSet&); - ~L1ECALPrefiringWeightProducer() override; - - static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); - -private: - void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override; - - double getPrefiringRate(double eta, double pt, TH2F* h_prefmap, fluctuations fluctuation) const; - - edm::EDGetTokenT > photons_token_; - edm::EDGetTokenT > jets_token_; - - TH2F* h_prefmap_photon; - TH2F* h_prefmap_jet; - std::string dataera_; - bool useEMpt_; - double prefiringRateSystUnc_; - double jetMaxMuonFraction_; - bool skipwarnings_; -}; - -L1ECALPrefiringWeightProducer::L1ECALPrefiringWeightProducer(const edm::ParameterSet& iConfig) { - photons_token_ = consumes >(iConfig.getParameter("ThePhotons")); - jets_token_ = consumes >(iConfig.getParameter("TheJets")); - - dataera_ = iConfig.getParameter("DataEra"); - useEMpt_ = iConfig.getParameter("UseJetEMPt"); - prefiringRateSystUnc_ = iConfig.getParameter("PrefiringRateSystematicUncty"); - jetMaxMuonFraction_ = iConfig.getParameter("JetMaxMuonFraction"); - skipwarnings_ = iConfig.getParameter("SkipWarnings"); - - TFile* file_prefiringmaps_; - std::string fname = iConfig.getParameter("L1Maps"); - edm::FileInPath mapsfilepath("PhysicsTools/PatUtils/data/" + fname); - file_prefiringmaps_ = new TFile(mapsfilepath.fullPath().c_str(), "read"); - if (file_prefiringmaps_ == nullptr && !skipwarnings_) - std::cout << "File with maps not found. All prefiring weights set to 0. " << std::endl; - - TString mapphotonfullname = "L1prefiring_photonptvseta_" + dataera_; - if (!file_prefiringmaps_->Get(mapphotonfullname) && !skipwarnings_) - std::cout << "Photon map not found. All photons prefiring weights set to 0. " << std::endl; - h_prefmap_photon = (TH2F*)file_prefiringmaps_->Get(mapphotonfullname); - - TString mapjetfullname = (useEMpt_) ? "L1prefiring_jetemptvseta_" + dataera_ : "L1prefiring_jetptvseta_" + dataera_; - if (!file_prefiringmaps_->Get(mapjetfullname) && !skipwarnings_) - std::cout << "Jet map not found. All jets prefiring weights set to 0. " << std::endl; - h_prefmap_jet = (TH2F*)file_prefiringmaps_->Get(mapjetfullname); - file_prefiringmaps_->Close(); - delete file_prefiringmaps_; - produces("nonPrefiringProb").setBranchAlias("nonPrefiringProb"); - produces("nonPrefiringProbUp").setBranchAlias("nonPrefiringProbUp"); - produces("nonPrefiringProbDown").setBranchAlias("nonPrefiringProbDown"); -} - -L1ECALPrefiringWeightProducer::~L1ECALPrefiringWeightProducer() { - delete h_prefmap_photon; - delete h_prefmap_jet; -} - -void L1ECALPrefiringWeightProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { - using namespace edm; - - //Photons - edm::Handle > thePhotons; - iEvent.getByToken(photons_token_, thePhotons); - - //Jets - edm::Handle > theJets; - iEvent.getByToken(jets_token_, theJets); - - //Probability for the event NOT to prefire, computed with the prefiring maps per object. - //Up and down values correspond to the resulting value when shifting up/down all prefiring rates in prefiring maps. - double nonPrefiringProba[3] = {1., 1., 1.}; //0: central, 1: up, 2: down - - for (const auto fluct : {fluctuations::central, fluctuations::up, fluctuations::down}) { - for (const auto& photon : *thePhotons) { - double pt_gam = photon.pt(); - double eta_gam = photon.eta(); - if (pt_gam < 20.) - continue; - if (fabs(eta_gam) < 2.) - continue; - if (fabs(eta_gam) > 3.) - continue; - double prefiringprob_gam = getPrefiringRate(eta_gam, pt_gam, h_prefmap_photon, fluct); - nonPrefiringProba[fluct] *= (1. - prefiringprob_gam); - } - - //Now applying the prefiring maps to jets in the affected regions. - for (const auto& jet : *theJets) { - double pt_jet = jet.pt(); - double eta_jet = jet.eta(); - double phi_jet = jet.phi(); - if (pt_jet < 20.) - continue; - if (fabs(eta_jet) < 2.) - continue; - if (fabs(eta_jet) > 3.) - continue; - if (jetMaxMuonFraction_ > 0 && jet.muonEnergyFraction() > jetMaxMuonFraction_) - continue; - //Loop over photons to remove overlap - double nonprefiringprobfromoverlappingphotons = 1.; - for (const auto& photon : *thePhotons) { - double pt_gam = photon.pt(); - double eta_gam = photon.eta(); - double phi_gam = photon.phi(); - if (pt_gam < 20.) - continue; - if (fabs(eta_gam) < 2.) - continue; - if (fabs(eta_gam) > 3.) - continue; - double dR = reco::deltaR(eta_jet, phi_jet, eta_gam, phi_gam); - if (dR > 0.4) - continue; - double prefiringprob_gam = getPrefiringRate(eta_gam, pt_gam, h_prefmap_photon, fluct); - nonprefiringprobfromoverlappingphotons *= (1. - prefiringprob_gam); - } - //useEMpt =true if one wants to use maps parametrized vs Jet EM pt instead of pt. - if (useEMpt_) - pt_jet *= (jet.neutralEmEnergyFraction() + jet.chargedEmEnergyFraction()); - double nonprefiringprobfromoverlappingjet = 1. - getPrefiringRate(eta_jet, pt_jet, h_prefmap_jet, fluct); - - if (nonprefiringprobfromoverlappingphotons == 1.) { - nonPrefiringProba[fluct] *= nonprefiringprobfromoverlappingjet; - } - //If overlapping photons have a non prefiring rate larger than the jet, then replace these weights by the jet one - else if (nonprefiringprobfromoverlappingphotons > nonprefiringprobfromoverlappingjet) { - if (nonprefiringprobfromoverlappingphotons != 0.) { - nonPrefiringProba[fluct] *= nonprefiringprobfromoverlappingjet / nonprefiringprobfromoverlappingphotons; - } else { - nonPrefiringProba[fluct] = 0.; - } - } - //Last case: if overlapping photons have a non prefiring rate smaller than the jet, don't consider the jet in the event weight, and do nothing. - } - } - - auto nonPrefiringProb = std::make_unique(nonPrefiringProba[0]); - auto nonPrefiringProbUp = std::make_unique(nonPrefiringProba[1]); - auto nonPrefiringProbDown = std::make_unique(nonPrefiringProba[2]); - iEvent.put(std::move(nonPrefiringProb), "nonPrefiringProb"); - iEvent.put(std::move(nonPrefiringProbUp), "nonPrefiringProbUp"); - iEvent.put(std::move(nonPrefiringProbDown), "nonPrefiringProbDown"); -} - -double L1ECALPrefiringWeightProducer::getPrefiringRate(double eta, - double pt, - TH2F* h_prefmap, - fluctuations fluctuation) const { - if (h_prefmap == nullptr && !skipwarnings_) - std::cout << "Prefiring map not found, setting prefiring rate to 0 " << std::endl; - if (h_prefmap == nullptr) - return 0.; - //Check pt is not above map overflow - int nbinsy = h_prefmap->GetNbinsY(); - double maxy = h_prefmap->GetYaxis()->GetBinLowEdge(nbinsy + 1); - if (pt >= maxy) - pt = maxy - 0.01; - int thebin = h_prefmap->FindBin(eta, pt); - - double prefrate = h_prefmap->GetBinContent(thebin); - - double statuncty = h_prefmap->GetBinError(thebin); - double systuncty = prefiringRateSystUnc_ * prefrate; - - if (fluctuation == up) - prefrate = std::min(1., prefrate + sqrt(pow(statuncty, 2) + pow(systuncty, 2))); - if (fluctuation == down) - prefrate = std::max(0., prefrate - sqrt(pow(statuncty, 2) + pow(systuncty, 2))); - return prefrate; -} - -// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ -void L1ECALPrefiringWeightProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { - edm::ParameterSetDescription desc; - - desc.add("ThePhotons", edm::InputTag("slimmedPhotons")); - desc.add("TheJets", edm::InputTag("slimmedJets")); - desc.add("L1Maps", "L1PrefiringMaps.root"); - desc.add("DataEra", "2017BtoF"); - desc.add("UseJetEMPt", false); - desc.add("PrefiringRateSystematicUncty", 0.2); - desc.add("JetMaxMuonFraction", 0.5); - desc.add("SkipWarnings", true); - descriptions.add("l1ECALPrefiringWeightProducer", desc); -} - -//define this as a plug-in -DEFINE_FWK_MODULE(L1ECALPrefiringWeightProducer); diff --git a/PhysicsTools/PatUtils/plugins/L1PrefiringWeightProducer.cc b/PhysicsTools/PatUtils/plugins/L1PrefiringWeightProducer.cc new file mode 100644 index 0000000000000..81850eb09cd50 --- /dev/null +++ b/PhysicsTools/PatUtils/plugins/L1PrefiringWeightProducer.cc @@ -0,0 +1,460 @@ +// -*- C++ -*- +// +// Package: ProdTutorial/L1PrefiringWeightProducer +// Class: L1PrefiringWeightProducer +// +/**\class L1PrefiringWeightProducer L1PrefiringWeightProducer.cc ProdTutorial/L1ECALPrefiringWeightProducer/plugins/L1ECALPrefiringWeightProducer.cc + + Description: [one line class summary] + + Implementation: + [Notes on implementation] +*/ +// +// Original Author: localusers user +// Created: Thu, 08 Nov 2018 16:16:00 GMT +// +// + +// system include files +#include + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" + +#include "DataFormats/PatCandidates/interface/Jet.h" +#include "DataFormats/PatCandidates/interface/Photon.h" +#include "DataFormats/PatCandidates/interface/Muon.h" + +#include "TFile.h" +#include "TF1.h" +#include "TH2.h" + +#include +enum fluctuations { central = 0, up, down, upStat, downStat, upSyst, downSyst }; + +class L1PrefiringWeightProducer : public edm::stream::EDProducer<> { +public: + explicit L1PrefiringWeightProducer(const edm::ParameterSet&); + ~L1PrefiringWeightProducer() override; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + void produce(edm::Event&, const edm::EventSetup&) override; + + double getPrefiringRateEcal(double eta, double pt, TH2F* h_prefmap, fluctuations fluctuation) const; + double getPrefiringRateMuon(double eta, double phi, double pt, fluctuations fluctuation) const; + + const edm::EDGetTokenT > photons_token_; + const edm::EDGetTokenT > jets_token_; + const edm::EDGetTokenT > muon_token_; + + const edm::EDPutTokenT nonPrefiringProbToken_; + const edm::EDPutTokenT nonPrefiringProbUpToken_; + const edm::EDPutTokenT nonPrefiringProbDownToken_; + + const edm::EDPutTokenT nonPrefiringProbECALToken_; + const edm::EDPutTokenT nonPrefiringProbECALUpToken_; + const edm::EDPutTokenT nonPrefiringProbECALDownToken_; + + const edm::EDPutTokenT nonPrefiringProbMuonToken_; + const edm::EDPutTokenT nonPrefiringProbMuonUpToken_; + const edm::EDPutTokenT nonPrefiringProbMuonDownToken_; + const edm::EDPutTokenT nonPrefiringProbMuonUpSystToken_; + const edm::EDPutTokenT nonPrefiringProbMuonDownSystToken_; + const edm::EDPutTokenT nonPrefiringProbMuonUpStatToken_; + const edm::EDPutTokenT nonPrefiringProbMuonDownStatToken_; + + std::unique_ptr file_prefiringmaps_; + std::unique_ptr file_prefiringparams_; + + TF1* parametrization0p0To0p2_; + TF1* parametrization0p2To0p3_; + TF1* parametrization0p3To0p55_; + TF1* parametrization0p55To0p83_; + TF1* parametrization0p83To1p24_; + TF1* parametrization1p24To1p4_; + TF1* parametrization1p4To1p6_; + TF1* parametrization1p6To1p8_; + TF1* parametrization1p8To2p1_; + TF1* parametrization2p1To2p25_; + TF1* parametrization2p25To2p4_; + TF1* parametrizationHotSpot_; + + TH2F* h_prefmap_photon_; + TH2F* h_prefmap_jet_; + const std::string dataeraEcal_; + const std::string dataeraMuon_; + const bool useEMpt_; + const double prefiringRateSystUncEcal_; + const double prefiringRateSystUncMuon_; + const double jetMaxMuonFraction_; + bool missingInputEcal_; + bool missingInputMuon_; +}; + +L1PrefiringWeightProducer::L1PrefiringWeightProducer(const edm::ParameterSet& iConfig) + : photons_token_(consumes >(iConfig.getParameter("ThePhotons"))), + jets_token_(consumes >(iConfig.getParameter("TheJets"))), + muon_token_(consumes >(iConfig.getParameter("TheMuons"))), + nonPrefiringProbToken_(produces("nonPrefiringProb")), + nonPrefiringProbUpToken_(produces("nonPrefiringProbUp")), + nonPrefiringProbDownToken_(produces("nonPrefiringProbDown")), + nonPrefiringProbECALToken_(produces("nonPrefiringProbECAL")), + nonPrefiringProbECALUpToken_(produces("nonPrefiringProbECALUp")), + nonPrefiringProbECALDownToken_(produces("nonPrefiringProbECALDown")), + nonPrefiringProbMuonToken_(produces("nonPrefiringProbMuon")), + nonPrefiringProbMuonUpToken_(produces("nonPrefiringProbMuonUp")), + nonPrefiringProbMuonDownToken_(produces("nonPrefiringProbMuonDown")), + nonPrefiringProbMuonUpSystToken_(produces("nonPrefiringProbMuonSystUp")), + nonPrefiringProbMuonDownSystToken_(produces("nonPrefiringProbMuonSystDown")), + nonPrefiringProbMuonUpStatToken_(produces("nonPrefiringProbMuonStatUp")), + nonPrefiringProbMuonDownStatToken_(produces("nonPrefiringProbMuonStatDown")), + dataeraEcal_(iConfig.getParameter("DataEraECAL")), + dataeraMuon_(iConfig.getParameter("DataEraMuon")), + useEMpt_(iConfig.getParameter("UseJetEMPt")), + prefiringRateSystUncEcal_(iConfig.getParameter("PrefiringRateSystematicUnctyECAL")), + prefiringRateSystUncMuon_(iConfig.getParameter("PrefiringRateSystematicUnctyMuon")), + jetMaxMuonFraction_(iConfig.getParameter("JetMaxMuonFraction")) { + missingInputEcal_ = false; + missingInputMuon_ = false; + + std::string fname = iConfig.getParameter("L1Maps"); + edm::FileInPath mapsfilepath("PhysicsTools/PatUtils/data/" + fname); + file_prefiringmaps_ = std::make_unique(mapsfilepath.fullPath().c_str(), "read"); + if (file_prefiringmaps_ == nullptr) { + missingInputEcal_ = true; + edm::LogError("L1PrefireWeightProducer") + << "File with maps not found. All prefiring weights set to 0. " << std::endl; + } + TString mapphotonfullname = "L1prefiring_photonptvseta_" + dataeraEcal_; + if (!file_prefiringmaps_->Get(mapphotonfullname)) { + missingInputEcal_ = true; + edm::LogError("L1PrefireWeightProducer") + << "Photon map not found. All photons prefiring weights set to 0. " << std::endl; + } + h_prefmap_photon_ = file_prefiringmaps_->Get(mapphotonfullname); + TString mapjetfullname = + (useEMpt_) ? "L1prefiring_jetemptvseta_" + dataeraEcal_ : "L1prefiring_jetptvseta_" + dataeraEcal_; + if (!file_prefiringmaps_->Get(mapjetfullname)) { + missingInputEcal_ = true; + edm::LogError("L1PrefireWeightProducer") << "Jet map not found. All jets prefiring weights set to 0. " << std::endl; + } + h_prefmap_jet_ = file_prefiringmaps_->Get(mapjetfullname); + file_prefiringmaps_->Close(); + + std::string fnameMuon = iConfig.getParameter("L1MuonParametrizations"); + edm::FileInPath paramsfilepath("PhysicsTools/PatUtils/data/" + fnameMuon); + file_prefiringparams_ = std::make_unique(paramsfilepath.fullPath().c_str(), "read"); + if (file_prefiringparams_ == nullptr) { + missingInputMuon_ = true; + edm::LogError("L1PrefireWeightProducer") + << "File with muon parametrizations not found. All prefiring weights set to 0." << std::endl; + } + TString paramName = "L1prefiring_muonparam_0.0To0.2_" + dataeraMuon_; + parametrization0p0To0p2_ = file_prefiringparams_->Get(paramName); + paramName = "L1prefiring_muonparam_0.2To0.3_" + dataeraMuon_; + parametrization0p2To0p3_ = file_prefiringparams_->Get(paramName); + paramName = "L1prefiring_muonparam_0.3To0.55_" + dataeraMuon_; + parametrization0p3To0p55_ = file_prefiringparams_->Get(paramName); + paramName = "L1prefiring_muonparam_0.55To0.83_" + dataeraMuon_; + parametrization0p55To0p83_ = file_prefiringparams_->Get(paramName); + paramName = "L1prefiring_muonparam_0.83To1.24_" + dataeraMuon_; + parametrization0p83To1p24_ = file_prefiringparams_->Get(paramName); + paramName = "L1prefiring_muonparam_1.24To1.4_" + dataeraMuon_; + parametrization1p24To1p4_ = file_prefiringparams_->Get(paramName); + paramName = "L1prefiring_muonparam_1.4To1.6_" + dataeraMuon_; + parametrization1p4To1p6_ = file_prefiringparams_->Get(paramName); + paramName = "L1prefiring_muonparam_1.6To1.8_" + dataeraMuon_; + parametrization1p6To1p8_ = file_prefiringparams_->Get(paramName); + paramName = "L1prefiring_muonparam_1.8To2.1_" + dataeraMuon_; + parametrization1p8To2p1_ = file_prefiringparams_->Get(paramName); + paramName = "L1prefiring_muonparam_2.1To2.25_" + dataeraMuon_; + parametrization2p1To2p25_ = file_prefiringparams_->Get(paramName); + paramName = "L1prefiring_muonparam_2.25To2.4_" + dataeraMuon_; + parametrization2p25To2p4_ = file_prefiringparams_->Get(paramName); + + if (parametrization0p0To0p2_ == nullptr || parametrization0p2To0p3_ == nullptr || + parametrization0p3To0p55_ == nullptr || parametrization0p55To0p83_ == nullptr || + parametrization0p83To1p24_ == nullptr || parametrization1p24To1p4_ == nullptr || + parametrization1p4To1p6_ == nullptr || parametrization1p6To1p8_ == nullptr || + parametrization1p8To2p1_ == nullptr || parametrization2p1To2p25_ == nullptr || + parametrization2p25To2p4_ == nullptr) { + missingInputMuon_ = true; + edm::LogError("L1PrefireWeightProducer") + << "Muon parametrization not found for at least one bin. All prefiring weights set to 0." << std::endl; + } + + paramName = "L1prefiring_muonparam_HotSpot_" + dataeraMuon_; + parametrizationHotSpot_ = file_prefiringparams_->Get(paramName); + file_prefiringparams_->Close(); + if ((dataeraMuon_.find("2016") != std::string::npos) && parametrizationHotSpot_ == nullptr) { + missingInputMuon_ = true; + edm::LogError("L1PrefireWeightProducer") + << "Year is 2016 and no Muon parametrization is found for hot spot. All prefiring weights set to 0." + << std::endl; + } +} + +L1PrefiringWeightProducer::~L1PrefiringWeightProducer() {} + +void L1PrefiringWeightProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + using namespace edm; + + //Photons + const std::vector& thePhotons = iEvent.get(photons_token_); + + //Jets + const std::vector& theJets = iEvent.get(jets_token_); + + //Muons + const std::vector& theMuons = iEvent.get(muon_token_); + + //Probability for the event NOT to prefire, computed with the prefiring maps per object. + //Up and down values correspond to the resulting value when shifting up/down all prefiring rates in prefiring maps. + double nonPrefiringProba[3] = {1., 1., 1.}; //0: central, 1: up, 2: down + double nonPrefiringProbaECAL[3] = {1., 1., 1.}; //0: central, 1: up, 2: down + double nonPrefiringProbaMuon[7] = { + 1., 1., 1., 1., 1., 1., 1.}; //0: central, 1: up, 2: down, 3: up stat, 4: down stat, 5: up syst, 6: down syst + + for (const auto fluct : {fluctuations::central, fluctuations::up, fluctuations::down}) { + if (!missingInputEcal_) { + for (const auto& photon : thePhotons) { + double pt_gam = photon.pt(); + double eta_gam = photon.eta(); + if (pt_gam < 20.) + continue; + if (fabs(eta_gam) < 2.) + continue; + if (fabs(eta_gam) > 3.) + continue; + double prefiringprob_gam = getPrefiringRateEcal(eta_gam, pt_gam, h_prefmap_photon_, fluct); + nonPrefiringProbaECAL[fluct] *= (1. - prefiringprob_gam); + } + + //Now applying the prefiring maps to jets in the affected regions. + for (const auto& jet : theJets) { + double pt_jet = jet.pt(); + double eta_jet = jet.eta(); + double phi_jet = jet.phi(); + if (pt_jet < 20.) + continue; + if (fabs(eta_jet) < 2.) + continue; + if (fabs(eta_jet) > 3.) + continue; + if (jetMaxMuonFraction_ > 0 && jet.muonEnergyFraction() > jetMaxMuonFraction_) + continue; + //Loop over photons to remove overlap + double nonprefiringprobfromoverlappingphotons = 1.; + bool foundOverlappingPhotons = false; + for (const auto& photon : thePhotons) { + double pt_gam = photon.pt(); + double eta_gam = photon.eta(); + double phi_gam = photon.phi(); + if (pt_gam < 20.) + continue; + if (fabs(eta_gam) < 2.) + continue; + if (fabs(eta_gam) > 3.) + continue; + double dR2 = reco::deltaR2(eta_jet, phi_jet, eta_gam, phi_gam); + if (dR2 > 0.16) + continue; + double prefiringprob_gam = getPrefiringRateEcal(eta_gam, pt_gam, h_prefmap_photon_, fluct); + nonprefiringprobfromoverlappingphotons *= (1. - prefiringprob_gam); + foundOverlappingPhotons = true; + } + //useEMpt =true if one wants to use maps parametrized vs Jet EM pt instead of pt. + if (useEMpt_) + pt_jet *= (jet.neutralEmEnergyFraction() + jet.chargedEmEnergyFraction()); + double nonprefiringprobfromoverlappingjet = 1. - getPrefiringRateEcal(eta_jet, pt_jet, h_prefmap_jet_, fluct); + + if (!foundOverlappingPhotons) { + nonPrefiringProbaECAL[fluct] *= nonprefiringprobfromoverlappingjet; + } + //If overlapping photons have a non prefiring rate larger than the jet, then replace these weights by the jet one + else if (nonprefiringprobfromoverlappingphotons > nonprefiringprobfromoverlappingjet) { + if (nonprefiringprobfromoverlappingphotons > 0.) { + nonPrefiringProbaECAL[fluct] *= nonprefiringprobfromoverlappingjet / nonprefiringprobfromoverlappingphotons; + } else { + nonPrefiringProbaECAL[fluct] = 0.; + } + } + //Last case: if overlapping photons have a non prefiring rate smaller than the jet, don't consider the jet in the event weight, and do nothing. + } + } + //Now calculate prefiring weights for muons + if (!missingInputMuon_) { + for (const auto& muon : theMuons) { + double pt = muon.pt(); + double phi = muon.phi(); + double eta = muon.eta(); + // Remove crappy tracker muons which would not have prefired the L1 trigger + if (pt < 5 || !muon.isLooseMuon()) + continue; + double prefiringprob_mu = getPrefiringRateMuon(eta, phi, pt, fluct); + nonPrefiringProbaMuon[fluct] *= (1. - prefiringprob_mu); + } + } + } + // Calculate combined weight as product of the weight for individual objects + for (const auto fluct : {fluctuations::central, fluctuations::up, fluctuations::down}) { + nonPrefiringProba[fluct] = nonPrefiringProbaECAL[fluct] * nonPrefiringProbaMuon[fluct]; + } + // Calculate statistical and systematic uncertainty separately in the muon case + for (const auto fluct : + {fluctuations::upSyst, fluctuations::downSyst, fluctuations::upStat, fluctuations::downStat}) { + if (!missingInputMuon_) { + for (const auto& muon : theMuons) { + double pt = muon.pt(); + double phi = muon.phi(); + double eta = muon.eta(); + // Remove crappy tracker muons which would not have prefired the L1 trigger + if (pt < 5 || !muon.isLooseMuon()) + continue; + double prefiringprob_mu = getPrefiringRateMuon(eta, phi, pt, fluct); + nonPrefiringProbaMuon[fluct] *= (1. - prefiringprob_mu); + } + } + } + //Move global prefire weights, as well as those for muons, photons, and jets, to the event + iEvent.emplace(nonPrefiringProbToken_, nonPrefiringProba[0]); + iEvent.emplace(nonPrefiringProbUpToken_, nonPrefiringProba[1]); + iEvent.emplace(nonPrefiringProbDownToken_, nonPrefiringProba[2]); + + iEvent.emplace(nonPrefiringProbECALToken_, nonPrefiringProbaECAL[0]); + iEvent.emplace(nonPrefiringProbECALUpToken_, nonPrefiringProbaECAL[1]); + iEvent.emplace(nonPrefiringProbECALDownToken_, nonPrefiringProbaECAL[2]); + + iEvent.emplace(nonPrefiringProbMuonToken_, nonPrefiringProbaMuon[0]); + iEvent.emplace(nonPrefiringProbMuonUpToken_, nonPrefiringProbaMuon[1]); + iEvent.emplace(nonPrefiringProbMuonDownToken_, nonPrefiringProbaMuon[2]); + iEvent.emplace(nonPrefiringProbMuonUpStatToken_, nonPrefiringProbaMuon[3]); + iEvent.emplace(nonPrefiringProbMuonDownStatToken_, nonPrefiringProbaMuon[4]); + iEvent.emplace(nonPrefiringProbMuonUpSystToken_, nonPrefiringProbaMuon[5]); + iEvent.emplace(nonPrefiringProbMuonDownSystToken_, nonPrefiringProbaMuon[6]); +} + +double L1PrefiringWeightProducer::getPrefiringRateEcal(double eta, + double pt, + TH2F* h_prefmap, + fluctuations fluctuation) const { + //Check pt is not above map overflow + int nbinsy = h_prefmap->GetNbinsY(); + double maxy = h_prefmap->GetYaxis()->GetBinLowEdge(nbinsy + 1); + if (pt >= maxy) + pt = maxy - 0.01; + int thebin = h_prefmap->FindBin(eta, pt); + + double prefrate = h_prefmap->GetBinContent(thebin); + + double statuncty = h_prefmap->GetBinError(thebin); + double systuncty = prefiringRateSystUncEcal_ * prefrate; + + if (fluctuation == up) + prefrate = std::min(1., prefrate + sqrt(pow(statuncty, 2) + pow(systuncty, 2))); + else if (fluctuation == down) + prefrate = std::max(0., prefrate - sqrt(pow(statuncty, 2) + pow(systuncty, 2))); + if (prefrate > 1.) { + edm::LogWarning("L1PrefireWeightProducer") << "Found a prefiring probability > 1. Setting to 1." << std::endl; + return 1.; + } + return prefrate; +} + +double L1PrefiringWeightProducer::getPrefiringRateMuon(double eta, + double phi, + double pt, + fluctuations fluctuation) const { + double prefrate; + double statuncty; + if ((dataeraMuon_.find("2016") != std::string::npos) && (eta > 1.24 && eta < 1.6) && + (phi > 2.44346 && phi < 2.79253)) { + prefrate = parametrizationHotSpot_->Eval(pt); + statuncty = parametrizationHotSpot_->GetParError(2); + } else if (std::abs(eta) < 0.2) { + prefrate = parametrization0p0To0p2_->Eval(pt); + statuncty = parametrization0p0To0p2_->GetParError(2); + } else if (std::abs(eta) < 0.3) { + prefrate = parametrization0p2To0p3_->Eval(pt); + statuncty = parametrization0p2To0p3_->GetParError(2); + } else if (std::abs(eta) < 0.55) { + prefrate = parametrization0p3To0p55_->Eval(pt); + statuncty = parametrization0p3To0p55_->GetParError(2); + } else if (std::abs(eta) < 0.83) { + prefrate = parametrization0p55To0p83_->Eval(pt); + statuncty = parametrization0p55To0p83_->GetParError(2); + } else if (std::abs(eta) < 1.24) { + prefrate = parametrization0p83To1p24_->Eval(pt); + statuncty = parametrization0p83To1p24_->GetParError(2); + } else if (std::abs(eta) < 1.4) { + prefrate = parametrization1p24To1p4_->Eval(pt); + statuncty = parametrization1p24To1p4_->GetParError(2); + } else if (std::abs(eta) < 1.6) { + prefrate = parametrization1p4To1p6_->Eval(pt); + statuncty = parametrization1p4To1p6_->GetParError(2); + } else if (std::abs(eta) < 1.8) { + prefrate = parametrization1p6To1p8_->Eval(pt); + statuncty = parametrization1p6To1p8_->GetParError(2); + } else if (std::abs(eta) < 2.1) { + prefrate = parametrization1p8To2p1_->Eval(pt); + statuncty = parametrization1p8To2p1_->GetParError(2); + } else if (std::abs(eta) < 2.25) { + prefrate = parametrization2p1To2p25_->Eval(pt); + statuncty = parametrization2p1To2p25_->GetParError(2); + } else if (std::abs(eta) < 2.4) { + prefrate = parametrization2p25To2p4_->Eval(pt); + statuncty = parametrization2p25To2p4_->GetParError(2); + } else { + LogDebug("L1PrefireWeightProducer") << "Muon outside of |eta| <= 2.4. Prefiring weight set to 0." << std::endl; + return 0.; + } + double systuncty = prefiringRateSystUncMuon_ * prefrate; + + if (fluctuation == up) + prefrate = std::min(1., prefrate + sqrt(pow(statuncty, 2) + pow(systuncty, 2))); + else if (fluctuation == down) + prefrate = std::max(0., prefrate - sqrt(pow(statuncty, 2) + pow(systuncty, 2))); + else if (fluctuation == upSyst) + prefrate = std::min(1., prefrate + systuncty); + else if (fluctuation == downSyst) + prefrate = std::max(0., prefrate - systuncty); + else if (fluctuation == upStat) + prefrate = std::min(1., prefrate + statuncty); + else if (fluctuation == downStat) + prefrate = std::max(0., prefrate - statuncty); + + if (prefrate > 1.) { + edm::LogWarning("L1PrefireWeightProducer") << "Found a prefiring probability > 1. Setting to 1." << std::endl; + return 1.; + } + return prefrate; +} + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void L1PrefiringWeightProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("TheMuons", edm::InputTag("slimmedMuons")); + desc.add("ThePhotons", edm::InputTag("slimmedPhotons")); + desc.add("TheJets", edm::InputTag("slimmedJets")); + desc.add("L1Maps", "L1PrefiringMaps.root"); + desc.add("L1MuonParametrizations", "L1MuonPrefiringParametriations.root"); + desc.add("DataEraECAL", "2017BtoF"); + desc.add("DataEraMuon", "2016"); + desc.add("UseJetEMPt", false); + desc.add("PrefiringRateSystematicUnctyECAL", 0.2); + desc.add("PrefiringRateSystematicUnctyMuon", 0.2); + desc.add("JetMaxMuonFraction", 0.5); + descriptions.add("l1PrefiringWeightProducer", desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(L1PrefiringWeightProducer); diff --git a/PhysicsTools/PatUtils/python/L1ECALPrefiringWeightProducer_cff.py b/PhysicsTools/PatUtils/python/L1ECALPrefiringWeightProducer_cff.py deleted file mode 100644 index c163719bd6fd4..0000000000000 --- a/PhysicsTools/PatUtils/python/L1ECALPrefiringWeightProducer_cff.py +++ /dev/null @@ -1,6 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -from PhysicsTools.PatUtils.l1ECALPrefiringWeightProducer_cfi import l1ECALPrefiringWeightProducer - -prefiringweight = l1ECALPrefiringWeightProducer.clone(SkipWarnings = False) - diff --git a/PhysicsTools/PatUtils/python/L1PrefiringWeightProducer_cff.py b/PhysicsTools/PatUtils/python/L1PrefiringWeightProducer_cff.py new file mode 100644 index 0000000000000..cb9ad47b5bcab --- /dev/null +++ b/PhysicsTools/PatUtils/python/L1PrefiringWeightProducer_cff.py @@ -0,0 +1,6 @@ +import FWCore.ParameterSet.Config as cms + +from PhysicsTools.PatUtils.l1PrefiringWeightProducer_cfi import l1PrefiringWeightProducer + +prefiringweight = l1PrefiringWeightProducer.clone() +