diff --git a/DataFormats/L1TCalorimeterPhase2/interface/CaloJet.h b/DataFormats/L1TCalorimeterPhase2/interface/CaloJet.h index d25ad2f686aca..5be863ec51179 100644 --- a/DataFormats/L1TCalorimeterPhase2/interface/CaloJet.h +++ b/DataFormats/L1TCalorimeterPhase2/interface/CaloJet.h @@ -20,7 +20,7 @@ namespace l1tp2 { inline float hovere() const { return hovere_; }; inline float isolation() const { return iso_; }; inline float puCorrPt() const { return puCorrPt_; }; - std::vector>& associated_l1EGs() { return associated_l1EGs_; }; + const std::vector>& associated_l1EGs() const { return associated_l1EGs_; }; void setExperimentalParams(const std::map& params) { experimentalParams_ = params; }; void setAssociated_l1EGs(const std::vector> l1EGs) { associated_l1EGs_ = l1EGs; }; diff --git a/L1Trigger/L1CaloTrigger/plugins/L1CaloJetHTTProducer.cc b/L1Trigger/L1CaloTrigger/plugins/L1CaloJetHTTProducer.cc new file mode 100644 index 0000000000000..d29032c46f497 --- /dev/null +++ b/L1Trigger/L1CaloTrigger/plugins/L1CaloJetHTTProducer.cc @@ -0,0 +1,123 @@ +// -*- C++ -*- +// +// Package: L1CaloTrigger +// Class: L1CaloJetHTTProducer +// +/**\class L1CaloJetHTTProducer L1CaloJetHTTProducer.cc + +Description: +Use the L1CaloJetProducer collections to calculate +HTT energy sum for CaloJets + +Implementation: +[Notes on implementation] +*/ +// +// Original Author: Tyler Ruggles +// Created: Fri Mar 22 2019 +// $Id$ +// +// + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include + +// Run2/PhaseI output formats +#include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h" +#include "DataFormats/L1Trigger/interface/Jet.h" +// GenJets if needed +#include "DataFormats/JetReco/interface/GenJet.h" +#include "DataFormats/JetReco/interface/GenJetCollection.h" + +class L1CaloJetHTTProducer : public edm::EDProducer { +public: + explicit L1CaloJetHTTProducer(const edm::ParameterSet&); + +private: + void produce(edm::Event&, const edm::EventSetup&) override; + + double EtaMax; + double PtMin; + + edm::EDGetTokenT> bxvCaloJetsToken_; + edm::Handle> bxvCaloJetsHandle; + + // Gen jet collections are only loaded and used if requested + // (use_gen_jets == true) + edm::EDGetTokenT> genJetsToken_; + edm::Handle> genJetsHandle; + + bool debug; + + bool use_gen_jets; +}; + +L1CaloJetHTTProducer::L1CaloJetHTTProducer(const edm::ParameterSet& iConfig) + : EtaMax(iConfig.getParameter("EtaMax")), + PtMin(iConfig.getParameter("PtMin")), + bxvCaloJetsToken_(consumes>(iConfig.getParameter("BXVCaloJetsInputTag"))), + genJetsToken_(consumes>(iConfig.getParameter("genJets"))), + debug(iConfig.getParameter("debug")), + use_gen_jets(iConfig.getParameter("use_gen_jets")) + +{ + produces("CaloJetHTT"); +} + +void L1CaloJetHTTProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + // Output collections + std::unique_ptr CaloJetHTT(new float); + + *CaloJetHTT = 0.; + + // CaloJet HTT for L1 collections + if (!use_gen_jets) { + iEvent.getByToken(bxvCaloJetsToken_, bxvCaloJetsHandle); + + if (bxvCaloJetsHandle.isValid()) { + for (const auto& caloJet : *bxvCaloJetsHandle.product()) { + if (caloJet.pt() < PtMin) + continue; + if (fabs(caloJet.eta()) > EtaMax) + continue; + *CaloJetHTT += float(caloJet.pt()); + } + } + + if (debug) { + LogDebug("L1CaloJetHTTProducer") << " BXV L1CaloJetCollection JetHTT = " << *CaloJetHTT << " for PtMin " << PtMin + << " and EtaMax " << EtaMax << "\n"; + } + } + + // CaloJet HTT for gen jets + if (use_gen_jets) { + iEvent.getByToken(genJetsToken_, genJetsHandle); + + if (genJetsHandle.isValid()) { + for (const auto& genJet : *genJetsHandle.product()) { + if (genJet.pt() < PtMin) + continue; + if (fabs(genJet.eta()) > EtaMax) + continue; + *CaloJetHTT += float(genJet.pt()); + } + } + + if (debug) { + LogDebug("L1CaloJetHTTProducer") << " Gen Jets HTT = " << *CaloJetHTT << " for PtMin " << PtMin << " and EtaMax " + << EtaMax << "\n"; + } + } + + iEvent.put(std::move(CaloJetHTT), "CaloJetHTT"); +} + +DEFINE_FWK_MODULE(L1CaloJetHTTProducer); diff --git a/L1Trigger/L1CaloTrigger/plugins/L1CaloJetProducer.cc b/L1Trigger/L1CaloTrigger/plugins/L1CaloJetProducer.cc new file mode 100644 index 0000000000000..9e99d6c99fb76 --- /dev/null +++ b/L1Trigger/L1CaloTrigger/plugins/L1CaloJetProducer.cc @@ -0,0 +1,1358 @@ +// -*- C++ -*- +// +// Package: L1CaloTrigger +// Class: L1CaloJetProducer +// +/**\class L1CaloJetProducer L1CaloJetProducer.cc + +Description: +Beginning with HCAL TPs, create HCAL jet, then +take L1EG crystal clusters from L1EGammaCrystalsProducer.cc +and clusters them within fixed number of trigger towers + +Implementation: +[Notes on implementation] +*/ +// +// Original Author: Tyler Ruggles +// Created: Tue Aug 29 2018 +// $Id$ +// +// + +#include "DataFormats/Math/interface/deltaR.h" +#include "DataFormats/Math/interface/deltaPhi.h" + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include + +#include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h" +//#include "DataFormats/L1TCalorimeterPhase2/interface/CaloCrystalCluster.h" +#include "DataFormats/L1TCalorimeterPhase2/interface/CaloJet.h" +#include "DataFormats/L1TCalorimeterPhase2/interface/CaloTower.h" +#include "DataFormats/L1THGCal/interface/HGCalTower.h" +#include "DataFormats/HcalDigi/interface/HcalDigiCollections.h" + +#include "CalibFormats/CaloTPG/interface/CaloTPGTranscoder.h" +#include "CalibFormats/CaloTPG/interface/CaloTPGRecord.h" + +#include "L1Trigger/L1TCalorimeter/interface/CaloTools.h" + +// For pT calibrations +#include "TF1.h" + +// Run2/PhaseI output formats +#include "DataFormats/L1Trigger/interface/Tau.h" +#include "DataFormats/L1Trigger/interface/Jet.h" + +class L1CaloJetProducer : public edm::EDProducer { +public: + explicit L1CaloJetProducer(const edm::ParameterSet &); + +private: + void produce(edm::Event &, const edm::EventSetup &) override; + //bool cluster_passes_base_cuts(float &cluster_pt, float &cluster_eta, float &iso, float &e2x5, float &e5x5) const; + int tower_diPhi(int &iPhi_1, int &iPhi_2) const; + int tower_diEta(int &iEta_1, int &iEta_2) const; + float get_deltaR(reco::Candidate::PolarLorentzVector &p4_1, reco::Candidate::PolarLorentzVector &p4_2) const; + float get_hcal_calibration(float &jet_pt, float &ecal_pt, float &ecal_L1EG_jet_pt, float &jet_eta) const; + float get_tau_pt_calibration(float &tau_pt, float &ecal_pt, float &l1EG_pt, float &n_L1EGs, float &tau_eta) const; + int loose_iso_tau_wp(float &tau_pt, float &tau_iso_et, float &tau_eta) const; + + double HcalTpEtMin; + double EcalTpEtMin; + double HGCalHadTpEtMin; + double HGCalEmTpEtMin; + double HFTpEtMin; + double EtMinForSeedHit; + double EtMinForCollection; + double EtMinForTauCollection; + + // For fetching jet calibrations + std::vector jetPtBins; + std::vector emFractionBinsBarrel; + std::vector absEtaBinsBarrel; + std::vector jetCalibrationsBarrel; + std::vector emFractionBinsHGCal; + std::vector absEtaBinsHGCal; + std::vector jetCalibrationsHGCal; + std::vector emFractionBinsHF; + std::vector absEtaBinsHF; + std::vector jetCalibrationsHF; + + // For fetching tau calibrations + std::vector tauPtBins; + std::vector tauAbsEtaBinsBarrel; + std::vector tauCalibrationsBarrel; + std::vector tauL1egInfoBarrel; + std::vector tauAbsEtaBinsHGCal; + std::vector tauCalibrationsHGCal; + std::vector tauL1egInfoHGCal; + + // For storing jet calibrations + std::vector>> calibrationsBarrel; + std::vector>> calibrationsHGCal; + std::vector>> calibrationsHF; + + // For storing tau calibration info + std::map> tauL1egInfoMapBarrel; + std::map> tauL1egInfoMapHGCal; + std::vector tauL1egValuesBarrel; // To preserve ordering + std::vector tauL1egValuesHGCal; // To preserve ordering + std::vector>>> tauPtCalibrationsBarrel; + std::vector>>> tauPtCalibrationsHGCal; + + bool debug; + edm::EDGetTokenT l1TowerToken_; + edm::Handle l1CaloTowerHandle; + + // TF1s defining tau isolation thresholds + TF1 isoTauBarrel = TF1("isoTauBarrelFunction", "([0] + [1]*TMath::Exp(-[2]*x))"); + TF1 isoTauHGCal = TF1("isoTauHGCalFunction", "([0] + [1]*TMath::Exp(-[2]*x))"); + + class l1CaloJetObj { + public: + bool barrelSeeded = true; // default to barrel seeded + reco::Candidate::PolarLorentzVector jetCluster; + reco::Candidate::PolarLorentzVector hcalJetCluster; + reco::Candidate::PolarLorentzVector ecalJetCluster; + reco::Candidate::PolarLorentzVector seedTower; + //reco::Candidate::PolarLorentzVector leadingL1EG; + reco::Candidate::PolarLorentzVector l1egJetCluster; + float jetClusterET = 0.; + float hcalJetClusterET = 0.; + float ecalJetClusterET = 0.; + float seedTowerET = 0.; + float leadingL1EGET = 0.; + float l1egJetClusterET = 0.; + + // For decay mode related checks with CaloTaus + std::vector> associated_l1EGs_; + + int seed_iEta = -99; + int seed_iPhi = -99; + + float hcal_seed = 0.; + //float hcal_3x3 = 0.; + float hcal_3x5 = 0.; + //float hcal_5x5 = 0.; + //float hcal_5x7 = 0.; + float hcal_7x7 = 0.; + //float hcal_2x3 = 0.; + //float hcal_2x3_1 = 0.; + //float hcal_2x3_2 = 0.; + //float hcal_2x2 = 0.; + //float hcal_2x2_1 = 0.; + //float hcal_2x2_2 = 0.; + //float hcal_2x2_3 = 0.; + //float hcal_2x2_4 = 0.; + float hcal_nHits = 0.; + + float ecal_seed = 0.; + //float ecal_3x3 = 0.; + float ecal_3x5 = 0.; + //float ecal_5x5 = 0.; + //float ecal_5x7 = 0.; + float ecal_7x7 = 0.; + //float ecal_2x3 = 0.; + //float ecal_2x3_1 = 0.; + //float ecal_2x3_2 = 0.; + //float ecal_2x2 = 0.; + //float ecal_2x2_1 = 0.; + //float ecal_2x2_2 = 0.; + //float ecal_2x2_3 = 0.; + //float ecal_2x2_4 = 0.; + float ecal_nHits = 0.; + + float l1eg_seed = 0.; + //float l1eg_3x3 = 0.; + float l1eg_3x5 = 0.; + //float l1eg_5x5 = 0.; + //float l1eg_5x7 = 0.; + float l1eg_7x7 = 0.; + //float l1eg_2x3 = 0.; + //float l1eg_2x3_1 = 0.; + //float l1eg_2x3_2 = 0.; + //float l1eg_2x2 = 0.; + //float l1eg_2x2_1 = 0.; + //float l1eg_2x2_2 = 0.; + //float l1eg_2x2_3 = 0.; + //float l1eg_2x2_4 = 0.; + float l1eg_nHits = 0.; + float n_l1eg_HoverE_LessThreshold = 0.; + + float l1eg_nL1EGs = 0.; + float l1eg_nL1EGs_standaloneSS = 0.; + float l1eg_nL1EGs_standaloneIso = 0.; + float l1eg_nL1EGs_trkMatchSS = 0.; + float l1eg_nL1EGs_trkMatchIso = 0.; + + float total_seed = 0.; + //float total_3x3 = 0.; + float total_3x5 = 0.; + //float total_5x5 = 0.; + //float total_5x7 = 0.; + float total_7x7 = 0.; + //float total_2x3 = 0.; + //float total_2x3_1 = 0.; + //float total_2x3_2 = 0.; + //float total_2x2 = 0.; + //float total_2x2_1 = 0.; + //float total_2x2_2 = 0.; + //float total_2x2_3 = 0.; + //float total_2x2_4 = 0.; + float total_nHits = 0.; + + void Init() { + SetJetClusterP4(0., 0., 0., 0.); + SetHcalJetClusterP4(0., 0., 0., 0.); + SetEcalJetClusterP4(0., 0., 0., 0.); + SetSeedP4(0., 0., 0., 0.); + //SetLeadingL1EGP4( 0., 0., 0., 0. ); + SetL1EGJetP4(0., 0., 0., 0.); + } + + void SetJetClusterP4(double pt, double eta, double phi, double mass) { + this->jetCluster.SetPt(pt); + this->jetCluster.SetEta(eta); + this->jetCluster.SetPhi(phi); + this->jetCluster.SetM(mass); + } + void SetHcalJetClusterP4(double pt, double eta, double phi, double mass) { + this->hcalJetCluster.SetPt(pt); + this->hcalJetCluster.SetEta(eta); + this->hcalJetCluster.SetPhi(phi); + this->hcalJetCluster.SetM(mass); + } + void SetEcalJetClusterP4(double pt, double eta, double phi, double mass) { + this->ecalJetCluster.SetPt(pt); + this->ecalJetCluster.SetEta(eta); + this->ecalJetCluster.SetPhi(phi); + this->ecalJetCluster.SetM(mass); + } + void SetSeedP4(double pt, double eta, double phi, double mass) { + this->seedTower.SetPt(pt); + this->seedTower.SetEta(eta); + this->seedTower.SetPhi(phi); + this->seedTower.SetM(mass); + } + //void SetLeadingL1EGP4( double pt, double eta, double phi, double mass ) + //{ + // this->leadingL1EG.SetPt( pt ); + // this->leadingL1EG.SetEta( eta ); + // this->leadingL1EG.SetPhi( phi ); + // this->leadingL1EG.SetM( mass ); + //} + void SetL1EGJetP4(double pt, double eta, double phi, double mass) { + this->l1egJetCluster.SetPt(pt); + this->l1egJetCluster.SetEta(eta); + this->l1egJetCluster.SetPhi(phi); + this->l1egJetCluster.SetM(mass); + } + }; + + class simpleL1obj { + public: + bool stale = false; // Hits become stale once used in clustering algorithm to prevent overlap in clusters + bool associated_with_tower = + false; // L1EGs become associated with a tower to find highest ET total for seeding jets + bool passesStandaloneSS = false; // Store whether any of the portions of a WP are passed + bool passesStandaloneIso = false; // Store whether any of the portions of a WP are passed + bool passesTrkMatchSS = false; // Store whether any of the portions of a WP are passed + bool passesTrkMatchIso = false; // Store whether any of the portions of a WP are passed + reco::Candidate::PolarLorentzVector p4; + + void SetP4(double pt, double eta, double phi, double mass) { + this->p4.SetPt(pt); + this->p4.SetEta(eta); + this->p4.SetPhi(phi); + this->p4.SetM(mass); + } + inline float pt() const { return p4.pt(); }; + inline float eta() const { return p4.eta(); }; + inline float phi() const { return p4.phi(); }; + inline float M() const { return p4.M(); }; + inline reco::Candidate::PolarLorentzVector GetP4() const { return p4; }; + }; + + class SimpleCaloHit { + public: + int towerIEta = -99; + int towerIPhi = -99; + float towerEta = -99; + float towerPhi = -99; + float ecalTowerEt = 0.; + float hcalTowerEt = 0.; + float l1egTowerEt = 0.; + float total_tower_et = 0.; + //float total_tower_plus_L1EGs_et=0.; + bool stale = false; // Hits become stale once used in clustering algorithm to prevent overlap in clusters + bool isBarrel = true; // Defaults to a barrel hit + + // L1EG info + int nL1eg = 0; + int l1egTrkSS = 0; + int l1egTrkIso = 0; + int l1egStandaloneSS = 0; + int l1egStandaloneIso = 0; + }; +}; + +L1CaloJetProducer::L1CaloJetProducer(const edm::ParameterSet &iConfig) + : HcalTpEtMin(iConfig.getParameter("HcalTpEtMin")), // Should default to 0 MeV + EcalTpEtMin(iConfig.getParameter("EcalTpEtMin")), // Should default to 0 MeV + HGCalHadTpEtMin(iConfig.getParameter("HGCalHadTpEtMin")), // Should default to 0 MeV + HGCalEmTpEtMin(iConfig.getParameter("HGCalEmTpEtMin")), // Should default to 0 MeV + HFTpEtMin(iConfig.getParameter("HFTpEtMin")), // Should default to 0 MeV + EtMinForSeedHit(iConfig.getParameter("EtMinForSeedHit")), // Should default to 2.5 GeV + EtMinForCollection(iConfig.getParameter("EtMinForCollection")), // Testing 10 GeV + EtMinForTauCollection(iConfig.getParameter("EtMinForTauCollection")), // Testing 10 GeV + jetPtBins(iConfig.getParameter>("jetPtBins")), + emFractionBinsBarrel(iConfig.getParameter>("emFractionBinsBarrel")), + absEtaBinsBarrel(iConfig.getParameter>("absEtaBinsBarrel")), + jetCalibrationsBarrel(iConfig.getParameter>("jetCalibrationsBarrel")), + emFractionBinsHGCal(iConfig.getParameter>("emFractionBinsHGCal")), + absEtaBinsHGCal(iConfig.getParameter>("absEtaBinsHGCal")), + jetCalibrationsHGCal(iConfig.getParameter>("jetCalibrationsHGCal")), + emFractionBinsHF(iConfig.getParameter>("emFractionBinsHF")), + absEtaBinsHF(iConfig.getParameter>("absEtaBinsHF")), + jetCalibrationsHF(iConfig.getParameter>("jetCalibrationsHF")), + tauPtBins(iConfig.getParameter>("tauPtBins")), + tauAbsEtaBinsBarrel(iConfig.getParameter>("tauAbsEtaBinsBarrel")), + tauCalibrationsBarrel(iConfig.getParameter>("tauCalibrationsBarrel")), + tauL1egInfoBarrel(iConfig.getParameter>("tauL1egInfoBarrel")), + tauAbsEtaBinsHGCal(iConfig.getParameter>("tauAbsEtaBinsHGCal")), + tauCalibrationsHGCal(iConfig.getParameter>("tauCalibrationsHGCal")), + tauL1egInfoHGCal(iConfig.getParameter>("tauL1egInfoHGCal")), + debug(iConfig.getParameter("debug")), + l1TowerToken_(consumes(iConfig.getParameter("l1CaloTowers"))) + +{ + produces("L1CaloJetsNoCuts"); + //produces("L1CaloJetsWithCuts"); + //produces("L1CaloClusterCollectionWithCuts"); + produces>("L1CaloJetCollectionBXV"); + produces>("L1CaloTauCollectionBXV"); + + if (debug) { + LogDebug("L1CaloJetProducer") << " HcalTpEtMin = " << HcalTpEtMin << " EcalTpEtMin = " << EcalTpEtMin << "\n"; + } + + // Fill the calibration 3D vector + // Dimension 1 is AbsEta bin + // Dimension 2 is EM Fraction bin + // Dimension 3 is jet pT bin which is filled with the actual callibration value + // size()-1 b/c the inputs have lower and upper bounds + // Do Barrel, then HGCal, then HF + int index = 0; + //calibrations[em_frac][abs_eta].push_back( jetCalibrationsBarrel.at(index) ); + for (unsigned int abs_eta = 0; abs_eta < absEtaBinsBarrel.size() - 1; abs_eta++) { + std::vector> em_bins; + for (unsigned int em_frac = 0; em_frac < emFractionBinsBarrel.size() - 1; em_frac++) { + std::vector pt_bin_calibs; + for (unsigned int pt = 0; pt < jetPtBins.size() - 1; pt++) { + pt_bin_calibs.push_back(jetCalibrationsBarrel.at(index)); + index++; + } + em_bins.push_back(pt_bin_calibs); + } + calibrationsBarrel.push_back(em_bins); + } + if (debug) { + LogDebug("L1CaloJetProducer") << " Loading Barrel calibrations: Loaded " << index + << " values vs. size() of input calibration file: " + << int(jetCalibrationsBarrel.size()) << "\n"; + } + + index = 0; + //calibrations[em_frac][abs_eta].push_back( jetCalibrationsHGCal.at(index) ); + for (unsigned int abs_eta = 0; abs_eta < absEtaBinsHGCal.size() - 1; abs_eta++) { + std::vector> em_bins; + for (unsigned int em_frac = 0; em_frac < emFractionBinsHGCal.size() - 1; em_frac++) { + std::vector pt_bin_calibs; + for (unsigned int pt = 0; pt < jetPtBins.size() - 1; pt++) { + pt_bin_calibs.push_back(jetCalibrationsHGCal.at(index)); + index++; + } + em_bins.push_back(pt_bin_calibs); + } + calibrationsHGCal.push_back(em_bins); + } + if (debug) { + LogDebug("L1CaloJetProducer") << " Loading HGCal calibrations: Loaded " << index + << " values vs. size() of input calibration file: " + << int(jetCalibrationsHGCal.size()) << "\n"; + } + + index = 0; + //calibrations[em_frac][abs_eta].push_back( jetCalibrationsHF.at(index) ); + for (unsigned int abs_eta = 0; abs_eta < absEtaBinsHF.size() - 1; abs_eta++) { + std::vector> em_bins; + for (unsigned int em_frac = 0; em_frac < emFractionBinsHF.size() - 1; em_frac++) { + std::vector pt_bin_calibs; + for (unsigned int pt = 0; pt < jetPtBins.size() - 1; pt++) { + pt_bin_calibs.push_back(jetCalibrationsHF.at(index)); + index++; + } + em_bins.push_back(pt_bin_calibs); + } + calibrationsHF.push_back(em_bins); + } + if (debug) { + LogDebug("L1CaloJetProducer") << " Loading HF calibrations: Loaded " << index + << " values vs. size() of input calibration file: " << int(jetCalibrationsHF.size()) + << "\n"; + } + + // Load Tau L1EG-base calibration info into maps + for (auto &first : tauL1egInfoBarrel) { + if (debug) { + LogDebug("L1CaloJetProducer") << " barrel l1egCount = " << first.getParameter("l1egCount") << "\n"; + for (auto &em_frac : first.getParameter>("l1egEmFractions")) { + LogDebug("L1CaloJetProducer") << " - EM = " << em_frac << "\n"; + } + } + double l1egCount = first.getParameter("l1egCount"); + std::vector l1egEmFractions = first.getParameter>("l1egEmFractions"); + tauL1egInfoMapBarrel[l1egCount] = l1egEmFractions; + tauL1egValuesBarrel.push_back(l1egCount); + } + std::sort(tauL1egValuesBarrel.begin(), tauL1egValuesBarrel.end()); + for (auto &first : tauL1egInfoHGCal) { + if (debug) { + LogDebug("L1CaloJetProducer") << " hgcal l1egCount = " << first.getParameter("l1egCount") << "\n"; + for (auto &em_frac : first.getParameter>("l1egEmFractions")) { + LogDebug("L1CaloJetProducer") << " - EM = " << em_frac << "\n"; + } + } + double l1egCount = first.getParameter("l1egCount"); + std::vector l1egEmFractions = first.getParameter>("l1egEmFractions"); + tauL1egInfoMapHGCal[l1egCount] = l1egEmFractions; + tauL1egValuesHGCal.push_back(l1egCount); + } + std::sort(tauL1egValuesHGCal.begin(), tauL1egValuesHGCal.end()); + // Fill the calibration 4D vector + // Dimension 1 is AbsEta bin + // Dimension 2 is L1EG count + // Dimension 3 is EM Fraction bin + // Dimension 4 is tau pT bin which is filled with the actual callibration value + // size()-1 b/c the inputs have lower and upper bounds (except L1EG b/c that is a cound) + // Do Barrel, then HGCal + // + // Note to future developers: be very concious of the order in which the calibrations are printed + // out in tool which makse the cfg files. You need to match that exactly when loading them and + // using the calibrations below. + index = 0; + for (unsigned int abs_eta = 0; abs_eta < tauAbsEtaBinsBarrel.size() - 1; abs_eta++) { + std::vector>> l1eg_bins; + for (auto &l1eg_info : tauL1egInfoMapBarrel) { + std::vector> em_bins; + for (unsigned int em_frac = 0; em_frac < l1eg_info.second.size() - 1; em_frac++) { + std::vector pt_bin_calibs; + for (unsigned int pt = 0; pt < tauPtBins.size() - 1; pt++) { + pt_bin_calibs.push_back(tauCalibrationsBarrel.at(index)); + index++; + } + em_bins.push_back(pt_bin_calibs); + } + l1eg_bins.push_back(em_bins); + } + tauPtCalibrationsBarrel.push_back(l1eg_bins); + } + if (debug) { + LogDebug("L1CaloJetProducer") << " Loading Barrel calibrations: Loaded " << index + << " values vs. size() of input calibration file: " + << int(tauCalibrationsBarrel.size()) << "\n"; + } + + index = 0; + for (unsigned int abs_eta = 0; abs_eta < tauAbsEtaBinsHGCal.size() - 1; abs_eta++) { + std::vector>> l1eg_bins; + for (const auto &l1eg_info : tauL1egInfoMapHGCal) { + std::vector> em_bins; + for (unsigned int em_frac = 0; em_frac < l1eg_info.second.size() - 1; em_frac++) { + std::vector pt_bin_calibs; + for (unsigned int pt = 0; pt < tauPtBins.size() - 1; pt++) { + pt_bin_calibs.push_back(tauCalibrationsHGCal.at(index)); + index++; + } + em_bins.push_back(pt_bin_calibs); + } + l1eg_bins.push_back(em_bins); + } + tauPtCalibrationsHGCal.push_back(l1eg_bins); + } + if (debug) { + LogDebug("L1CaloJetProducer") << " Loading HGCal calibrations: Loaded " << index + << " values vs. size() of input calibration file: " + << int(tauCalibrationsHGCal.size()) << "\n"; + } + + isoTauBarrel.SetParameter(0, 0.25); + isoTauBarrel.SetParameter(1, 0.85); + isoTauBarrel.SetParameter(2, 0.094); + isoTauHGCal.SetParameter(0, 0.34); + isoTauHGCal.SetParameter(1, 0.35); + isoTauHGCal.SetParameter(2, 0.051); +} + +void L1CaloJetProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) { + // Output collections + std::unique_ptr L1CaloJetsNoCuts(new l1tp2::CaloJetsCollection); + //std::unique_ptr L1CaloJetsWithCuts( new l1tp2::CaloJetsCollection ); + //std::unique_ptr L1CaloClusterCollectionWithCuts( new l1extra::L1JetParticleCollection ); + std::unique_ptr> L1CaloJetCollectionBXV(new l1t::JetBxCollection); + std::unique_ptr> L1CaloTauCollectionBXV(new l1t::TauBxCollection); + + // Load the ECAL+HCAL tower sums coming from L1EGammaCrystalsEmulatorProducer.cc + std::vector l1CaloTowers; + + iEvent.getByToken(l1TowerToken_, l1CaloTowerHandle); + for (auto &hit : *l1CaloTowerHandle.product()) { + SimpleCaloHit l1Hit; + l1Hit.ecalTowerEt = hit.ecalTowerEt(); + l1Hit.hcalTowerEt = hit.hcalTowerEt(); + l1Hit.l1egTowerEt = hit.l1egTowerEt(); + // Add min ET thresholds for tower ET + if (l1Hit.ecalTowerEt < EcalTpEtMin) + l1Hit.ecalTowerEt = 0.0; + if (l1Hit.hcalTowerEt < HcalTpEtMin) + l1Hit.hcalTowerEt = 0.0; + l1Hit.total_tower_et = l1Hit.ecalTowerEt + l1Hit.hcalTowerEt + l1Hit.l1egTowerEt; + l1Hit.towerIEta = hit.towerIEta(); + l1Hit.towerIPhi = hit.towerIPhi(); + l1Hit.nL1eg = hit.nL1eg(); + l1Hit.l1egTrkSS = hit.l1egTrkSS(); + l1Hit.l1egTrkIso = hit.l1egTrkIso(); + l1Hit.l1egStandaloneSS = hit.l1egStandaloneSS(); + l1Hit.l1egStandaloneIso = hit.l1egStandaloneIso(); + l1Hit.isBarrel = hit.isBarrel(); + + // FIXME There is an error in the L1EGammaCrystalsEmulatorProducer.cc which is + // returning towers with minimal ECAL energy, and no HCAL energy with these + // iEta/iPhi coordinates and eta = -88.653152 and phi = -99.000000. + // Skip these for the time being until the upstream code has been debugged + if ((int)l1Hit.towerIEta == -1016 && (int)l1Hit.towerIPhi == -962) + continue; + + l1Hit.towerEta = hit.towerEta(); + l1Hit.towerPhi = hit.towerPhi(); + l1CaloTowers.push_back(l1Hit); + if (debug) { + LogDebug("L1CaloJetProducer") << " Tower iEta " << (int)l1Hit.towerIEta << " iPhi " << (int)l1Hit.towerIPhi + << " eta " << l1Hit.towerEta << " phi " << l1Hit.towerPhi << " ecal_et " + << l1Hit.ecalTowerEt << " hacl_et " << l1Hit.hcalTowerEt << " total_et " + << l1Hit.total_tower_et << "\n"; + } + } + + // Sort the ECAL+HCAL+L1EGs tower sums based on total ET + std::sort(begin(l1CaloTowers), end(l1CaloTowers), [](const SimpleCaloHit &a, SimpleCaloHit &b) { + return a.total_tower_et > b.total_tower_et; + }); + + /************************************************************************** + * Begin with making CaloJets in 9x9 grid based on all energy not included in L1EG Objs. + * For reference, Run-I used 12x12 grid and Stage-2/Phase-I used 9x9 grid. + * We plan to further study this choice and possibly move towards a more circular shape + * Create jetCluster within 9x9 of highest ET seed tower. + * 9 trigger towers contains all of an ak-0.4 jets, but overshoots on the corners. + ******************************************************************************/ + + // Experimental parameters, don't want to bother with hardcoding them in data format + std::map params; + + std::vector l1CaloJetObjs; + + // Count the number of unused HCAL TPs so we can stop while loop after done. + // Clustering can also stop once there are no seed hits >= EtMinForSeedHit + int n_towers = l1CaloTowers.size(); + int n_stale = 0; + bool caloJetClusteringFinished = false; + while (!caloJetClusteringFinished && n_towers != n_stale) { + l1CaloJetObj caloJetObj; + caloJetObj.Init(); + + // First find highest ET ECAL+HCAL+L1EGs tower and use to seed the 9x9 Jet + int cnt = 0; + for (auto &l1CaloTower : l1CaloTowers) { + cnt++; + if (l1CaloTower.stale) + continue; // skip l1CaloTowers which are already used + + if (caloJetObj.jetClusterET == 0.0) // this is the first l1CaloTower to seed the jet + { + // Check if the leading unused tower has ET < min for seeding a jet. + // If so, stop jet clustering + if (l1CaloTower.total_tower_et < EtMinForSeedHit) { + caloJetClusteringFinished = true; + continue; + } + l1CaloTower.stale = true; + n_stale++; + + // Set seed location needed for delta iEta/iPhi, eta/phi comparisons later + if (l1CaloTower.isBarrel) + caloJetObj.barrelSeeded = true; + else + caloJetObj.barrelSeeded = false; + + // 3 4-vectors for ECAL, HCAL, ECAL+HCAL for adding together + reco::Candidate::PolarLorentzVector hcalP4( + l1CaloTower.hcalTowerEt, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.); + reco::Candidate::PolarLorentzVector ecalP4( + l1CaloTower.ecalTowerEt, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.); + reco::Candidate::PolarLorentzVector l1egP4( + l1CaloTower.l1egTowerEt, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.); + reco::Candidate::PolarLorentzVector totalP4( + l1CaloTower.total_tower_et, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.); + + if (hcalP4.pt() > 0) { + caloJetObj.hcal_nHits++; + caloJetObj.hcalJetCluster += hcalP4; + caloJetObj.hcalJetClusterET += l1CaloTower.hcalTowerEt; + } + if (ecalP4.pt() > 0) { + caloJetObj.ecal_nHits++; + caloJetObj.ecalJetCluster += ecalP4; + caloJetObj.ecalJetClusterET += l1CaloTower.ecalTowerEt; + } + if (l1egP4.pt() > 0) { + caloJetObj.l1eg_nHits++; + caloJetObj.l1egJetCluster += l1egP4; + caloJetObj.l1egJetClusterET += l1CaloTower.l1egTowerEt; + caloJetObj.l1eg_nL1EGs += l1CaloTower.nL1eg; + + caloJetObj.l1eg_nL1EGs_standaloneSS += l1CaloTower.l1egStandaloneSS; + caloJetObj.l1eg_nL1EGs_standaloneIso += l1CaloTower.l1egStandaloneIso; + caloJetObj.l1eg_nL1EGs_trkMatchSS += l1CaloTower.l1egTrkSS; + caloJetObj.l1eg_nL1EGs_trkMatchIso += l1CaloTower.l1egTrkIso; + + if (l1CaloTower.isBarrel) { + // For decay mode related checks with CaloTaus + // only applicable in the barrel at the moment: + // l1eg pt, HCAL ET, ECAL ET, dEta, dPhi, trkSS, trkIso, standaloneSS, standaloneIso + std::vector l1EG_info = {float(l1egP4.pt()), + float(hcalP4.pt()), + float(ecalP4.pt()), + 0., + 0., + float(l1CaloTower.l1egTrkSS), + float(l1CaloTower.l1egTrkIso), + float(l1CaloTower.l1egStandaloneSS), + float(l1CaloTower.l1egStandaloneIso)}; + if (l1EG_info[1] / (l1EG_info[0] + l1EG_info[2]) < 0.25) { + caloJetObj.n_l1eg_HoverE_LessThreshold++; + } + caloJetObj.associated_l1EGs_.push_back(l1EG_info); + } + } + if (totalP4.pt() > 0) { + caloJetObj.total_nHits++; + caloJetObj.jetCluster += totalP4; + caloJetObj.jetClusterET += l1CaloTower.total_tower_et; + caloJetObj.seedTower += totalP4; + caloJetObj.seedTowerET += l1CaloTower.total_tower_et; + } + + caloJetObj.seed_iEta = l1CaloTower.towerIEta; + caloJetObj.seed_iPhi = l1CaloTower.towerIPhi; + + if (debug) { + LogDebug("L1CaloJetProducer") << " -- hit " << cnt << " , seeding input p4 pt " + << l1CaloTower.total_tower_et << " eta " << l1CaloTower.towerEta << " phi " + << l1CaloTower.towerPhi << "\n"; + LogDebug("L1CaloJetProducer") << " -- hit " << cnt << " , seeding input2 p4 pt " << totalP4.pt() << " eta " + << totalP4.eta() << " phi " << totalP4.phi() << "\n"; + LogDebug("L1CaloJetProducer") << " -- hit " << cnt << " seeding reslt tot p4 pt " << caloJetObj.jetClusterET + << " eta " << caloJetObj.jetCluster.eta() << " phi " + << caloJetObj.jetCluster.phi() << "\n"; + } + + // Need to add the seed energy to the dR rings + caloJetObj.hcal_seed += hcalP4.pt(); + //caloJetObj.hcal_3x3 += hcalP4.pt(); + caloJetObj.hcal_3x5 += hcalP4.pt(); + //caloJetObj.hcal_5x5 += hcalP4.pt(); + //caloJetObj.hcal_5x7 += hcalP4.pt(); + caloJetObj.hcal_7x7 += hcalP4.pt(); + caloJetObj.ecal_seed += ecalP4.pt(); + //caloJetObj.ecal_3x3 += ecalP4.pt(); + caloJetObj.ecal_3x5 += ecalP4.pt(); + //caloJetObj.ecal_5x5 += ecalP4.pt(); + //caloJetObj.ecal_5x7 += ecalP4.pt(); + caloJetObj.ecal_7x7 += ecalP4.pt(); + caloJetObj.l1eg_seed += l1egP4.pt(); + //caloJetObj.l1eg_3x3 += l1egP4.pt(); + caloJetObj.l1eg_3x5 += l1egP4.pt(); + //caloJetObj.l1eg_5x5 += l1egP4.pt(); + //caloJetObj.l1eg_5x7 += l1egP4.pt(); + caloJetObj.l1eg_7x7 += l1egP4.pt(); + caloJetObj.total_seed += totalP4.pt(); + //caloJetObj.total_3x3 += totalP4.pt(); + caloJetObj.total_3x5 += totalP4.pt(); + //caloJetObj.total_5x5 += totalP4.pt(); + //caloJetObj.total_5x7 += totalP4.pt(); + caloJetObj.total_7x7 += totalP4.pt(); + + // Some discrimination vars, 2x2s and 2x3 including central seed + //caloJetObj.hcal_2x3 += hcalP4.pt(); + //caloJetObj.hcal_2x3_1 += hcalP4.pt(); + //caloJetObj.hcal_2x3_2 += hcalP4.pt(); + //caloJetObj.ecal_2x3 += ecalP4.pt(); + //caloJetObj.ecal_2x3_1 += ecalP4.pt(); + //caloJetObj.ecal_2x3_2 += ecalP4.pt(); + //caloJetObj.l1eg_2x3 += l1egP4.pt(); + //caloJetObj.l1eg_2x3_1 += l1egP4.pt(); + //caloJetObj.l1eg_2x3_2 += l1egP4.pt(); + //caloJetObj.total_2x3 += totalP4.pt(); + //caloJetObj.total_2x3_1 += totalP4.pt(); + //caloJetObj.total_2x3_2 += totalP4.pt(); + + //caloJetObj.hcal_2x2 += hcalP4.pt(); + //caloJetObj.hcal_2x2_1 += hcalP4.pt(); + //caloJetObj.hcal_2x2_2 += hcalP4.pt(); + //caloJetObj.hcal_2x2_3 += hcalP4.pt(); + //caloJetObj.hcal_2x2_4 += hcalP4.pt(); + //caloJetObj.ecal_2x2 += ecalP4.pt(); + //caloJetObj.ecal_2x2_1 += ecalP4.pt(); + //caloJetObj.ecal_2x2_2 += ecalP4.pt(); + //caloJetObj.ecal_2x2_3 += ecalP4.pt(); + //caloJetObj.ecal_2x2_4 += ecalP4.pt(); + //caloJetObj.l1eg_2x2 += l1egP4.pt(); + //caloJetObj.l1eg_2x2_1 += l1egP4.pt(); + //caloJetObj.l1eg_2x2_2 += l1egP4.pt(); + //caloJetObj.l1eg_2x2_3 += l1egP4.pt(); + //caloJetObj.l1eg_2x2_4 += l1egP4.pt(); + //caloJetObj.total_2x2 += totalP4.pt(); + //caloJetObj.total_2x2_1 += totalP4.pt(); + //caloJetObj.total_2x2_2 += totalP4.pt(); + //caloJetObj.total_2x2_3 += totalP4.pt(); + //caloJetObj.total_2x2_4 += totalP4.pt(); + continue; + } + + // Unused l1CaloTowers which are not the initial seed + // Depending on seed and tower locations calculate iEta/iPhi or eta/phi comparisons. + // The defaults of 99 will automatically fail comparisons for the incorrect regions. + int hit_iPhi = 99; + int d_iEta = 99; + int d_iPhi = 99; + float d_eta = 99; + float d_phi = 99; + if (caloJetObj.barrelSeeded && l1CaloTower.isBarrel) // use iEta/iPhi comparisons + { + hit_iPhi = l1CaloTower.towerIPhi; + d_iEta = tower_diEta(caloJetObj.seed_iEta, l1CaloTower.towerIEta); + d_iPhi = tower_diPhi(caloJetObj.seed_iPhi, hit_iPhi); + } else // either seed or tower are in HGCal or HF, use eta/phi + { + d_eta = caloJetObj.seedTower.eta() - l1CaloTower.towerEta; + d_phi = reco::deltaPhi(caloJetObj.seedTower.phi(), l1CaloTower.towerPhi); + } + + // 7x7 HCAL Trigger Towers + // If seeded in barrel and hit is barrel then we can compare iEta/iPhi, else need to use eta/phi + // in HGCal / transition region + if ((abs(d_iEta) <= 3 && abs(d_iPhi) <= 3) || (fabs(d_eta) < 0.3 && fabs(d_phi) < 0.3)) { + l1CaloTower.stale = true; + n_stale++; + + // 3 4-vectors for ECAL, HCAL, ECAL+HCAL for adding together + reco::Candidate::PolarLorentzVector hcalP4( + l1CaloTower.hcalTowerEt, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.); + reco::Candidate::PolarLorentzVector ecalP4( + l1CaloTower.ecalTowerEt, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.); + reco::Candidate::PolarLorentzVector l1egP4( + l1CaloTower.l1egTowerEt, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.); + reco::Candidate::PolarLorentzVector totalP4( + l1CaloTower.total_tower_et, l1CaloTower.towerEta, l1CaloTower.towerPhi, 0.); + + if (hcalP4.pt() > 0) { + caloJetObj.hcal_nHits++; + caloJetObj.hcalJetCluster += hcalP4; + caloJetObj.hcalJetClusterET += l1CaloTower.hcalTowerEt; + } + if (ecalP4.pt() > 0) { + caloJetObj.ecal_nHits++; + caloJetObj.ecalJetCluster += ecalP4; + caloJetObj.ecalJetClusterET += l1CaloTower.ecalTowerEt; + } + if (l1egP4.pt() > 0) { + caloJetObj.l1eg_nHits++; + caloJetObj.l1egJetCluster += l1egP4; + caloJetObj.l1egJetClusterET += l1CaloTower.l1egTowerEt; + caloJetObj.l1eg_nL1EGs += l1CaloTower.nL1eg; + } + if (totalP4.pt() > 0) { + caloJetObj.total_nHits++; + caloJetObj.jetCluster += totalP4; + caloJetObj.jetClusterET += l1CaloTower.total_tower_et; + } + + if (debug) { + LogDebug("L1CaloJetProducer") << " ---- hit " << cnt << " input p4 pt " << totalP4.pt() << " eta " + << totalP4.eta() << " phi " << totalP4.phi() << "\n"; + LogDebug("L1CaloJetProducer") << " ---- hit " << cnt << " resulting p4 pt " << caloJetObj.jetClusterET + << " eta " << caloJetObj.jetCluster.eta() << " phi " + << caloJetObj.jetCluster.phi() << "\n"; + } + + if ((abs(d_iEta) == 0 && abs(d_iPhi) == 0) || (fabs(d_eta) < 0.043 && fabs(d_phi) < 0.043)) { + caloJetObj.hcal_seed += hcalP4.pt(); + caloJetObj.ecal_seed += ecalP4.pt(); + caloJetObj.l1eg_seed += l1egP4.pt(); + caloJetObj.total_seed += totalP4.pt(); + } + //if ( (abs( d_iEta ) <= 1 && abs( d_iPhi ) <= 1) || + // ( fabs( d_eta ) < 0.13 && fabs( d_phi ) < 0.13 ) ) + //{ + // caloJetObj.hcal_3x3 += hcalP4.pt(); + // caloJetObj.ecal_3x3 += ecalP4.pt(); + // caloJetObj.l1eg_3x3 += l1egP4.pt(); + // caloJetObj.total_3x3 += totalP4.pt(); + //} + if ((abs(d_iEta) <= 1 && abs(d_iPhi) <= 2) || (fabs(d_eta) < 0.13 && fabs(d_phi) < 0.22)) { + caloJetObj.hcal_3x5 += hcalP4.pt(); + caloJetObj.ecal_3x5 += ecalP4.pt(); + caloJetObj.l1eg_3x5 += l1egP4.pt(); + caloJetObj.total_3x5 += totalP4.pt(); + + // Do this for 3x5 only + if (l1egP4.pt() > 0) { + caloJetObj.l1eg_nL1EGs_standaloneSS += l1CaloTower.l1egStandaloneSS; + caloJetObj.l1eg_nL1EGs_standaloneIso += l1CaloTower.l1egStandaloneIso; + caloJetObj.l1eg_nL1EGs_trkMatchSS += l1CaloTower.l1egTrkSS; + caloJetObj.l1eg_nL1EGs_trkMatchIso += l1CaloTower.l1egTrkIso; + + // For decay mode related checks with CaloTaus + // only applicable in the barrel at the moment: + // l1eg pt, HCAL ET, ECAL ET, d_iEta, d_iPhi, trkSS, trkIso, standaloneSS, standaloneIso + std::vector l1EG_info = {float(l1egP4.pt()), + float(hcalP4.pt()), + float(ecalP4.pt()), + float(d_iEta), + float(d_iPhi), + float(l1CaloTower.l1egTrkSS), + float(l1CaloTower.l1egTrkIso), + float(l1CaloTower.l1egStandaloneSS), + float(l1CaloTower.l1egStandaloneIso)}; + if (l1EG_info[1] / (l1EG_info[0] + l1EG_info[2]) < 0.25) { + caloJetObj.n_l1eg_HoverE_LessThreshold++; + } + caloJetObj.associated_l1EGs_.push_back(l1EG_info); + } + } + //if ( ( abs( d_iEta ) <= 2 && abs( d_iPhi ) <= 2) || + // ( fabs( d_eta ) < 0.22 && fabs( d_phi ) < 0.22 ) ) + //{ + // caloJetObj.hcal_5x5 += hcalP4.pt(); + // caloJetObj.ecal_5x5 += ecalP4.pt(); + // caloJetObj.l1eg_5x5 += l1egP4.pt(); + // caloJetObj.total_5x5 += totalP4.pt(); + //} + //if ( ( abs( d_iEta ) <= 2 && abs( d_iPhi ) <= 3) || + // ( fabs( d_eta ) < 0.22 && fabs( d_phi ) < 0.3 ) ) + //{ + // caloJetObj.hcal_5x7 += hcalP4.pt(); + // caloJetObj.ecal_5x7 += ecalP4.pt(); + // caloJetObj.l1eg_5x7 += l1egP4.pt(); + // caloJetObj.total_5x7 += totalP4.pt(); + //} + if ((abs(d_iEta) <= 3 && abs(d_iPhi) <= 3) || (fabs(d_eta) < 0.3 && fabs(d_phi) < 0.3)) { + caloJetObj.hcal_7x7 += hcalP4.pt(); + caloJetObj.ecal_7x7 += ecalP4.pt(); + caloJetObj.l1eg_7x7 += l1egP4.pt(); + caloJetObj.total_7x7 += totalP4.pt(); + } + + //// Some discrimination vars, 2x2s and 2x3s including central seed + //// Barrel first, 2x3 + //if ( ( d_iEta == 0 || d_iEta == 1 ) && abs(d_iPhi) <= 1 ) + //{ + // caloJetObj.hcal_2x3_1 += hcalP4.pt(); + // caloJetObj.ecal_2x3_1 += ecalP4.pt(); + // caloJetObj.l1eg_2x3_1 += l1egP4.pt(); + // caloJetObj.total_2x3_1 += totalP4.pt(); + //} + //if ( ( d_iEta == 0 || d_iEta == -1 ) && abs(d_iPhi) <= 1 ) + //{ + // caloJetObj.hcal_2x3_2 += hcalP4.pt(); + // caloJetObj.ecal_2x3_2 += ecalP4.pt(); + // caloJetObj.l1eg_2x3_2 += l1egP4.pt(); + // caloJetObj.total_2x3_2 += totalP4.pt(); + //} + //// HGCal / HF + //if ( fabs( d_eta ) < 0.087 && fabs( d_phi ) < 0.13 ) + //{ + // caloJetObj.hcal_2x3 += hcalP4.pt(); + // caloJetObj.ecal_2x3 += ecalP4.pt(); + // caloJetObj.l1eg_2x3 += l1egP4.pt(); + // caloJetObj.total_2x3 += totalP4.pt(); + //} + + //// Now 2x2 Barrel first + //if ( ( d_iEta == 0 || d_iEta == 1 ) && ( d_iPhi == 0 || d_iPhi == 1 ) ) + //{ + // caloJetObj.hcal_2x2_1 += hcalP4.pt(); + // caloJetObj.ecal_2x2_1 += ecalP4.pt(); + // caloJetObj.l1eg_2x2_1 += l1egP4.pt(); + // caloJetObj.total_2x2_1 += totalP4.pt(); + //} + //if ( ( d_iEta == 0 || d_iEta == 1 ) && ( d_iPhi == 0 || d_iPhi == -1 ) ) + //{ + // caloJetObj.hcal_2x2_2 += hcalP4.pt(); + // caloJetObj.ecal_2x2_2 += ecalP4.pt(); + // caloJetObj.l1eg_2x2_2 += l1egP4.pt(); + // caloJetObj.total_2x2_2 += totalP4.pt(); + //} + //if ( ( d_iEta == 0 || d_iEta == -1 ) && ( d_iPhi == 0 || d_iPhi == 1 ) ) + //{ + // caloJetObj.hcal_2x2_3 += hcalP4.pt(); + // caloJetObj.ecal_2x2_3 += ecalP4.pt(); + // caloJetObj.l1eg_2x2_3 += l1egP4.pt(); + // caloJetObj.total_2x2_3 += totalP4.pt(); + //} + //if ( ( d_iEta == 0 || d_iEta == -1 ) && ( d_iPhi == 0 || d_iPhi == -1 ) ) + //{ + // caloJetObj.hcal_2x2_4 += hcalP4.pt(); + // caloJetObj.ecal_2x2_4 += ecalP4.pt(); + // caloJetObj.l1eg_2x2_4 += l1egP4.pt(); + // caloJetObj.total_2x2_4 += totalP4.pt(); + //} + //// HGCal / HF + //if ( fabs( d_eta ) < 0.087 && fabs( d_phi ) < 0.087 ) + //{ + // caloJetObj.hcal_2x2 += hcalP4.pt(); + // caloJetObj.ecal_2x2 += ecalP4.pt(); + // caloJetObj.l1eg_2x2 += l1egP4.pt(); + // caloJetObj.total_2x2 += totalP4.pt(); + //} + } + } + + if (caloJetObj.jetClusterET > 0.0) { + l1CaloJetObjs.push_back(caloJetObj); + } + + } // end while loop of HCAL TP clustering + + // Sort JetClusters so we can begin with the highest pt for next step of jet clustering + std::sort(begin(l1CaloJetObjs), end(l1CaloJetObjs), [](const l1CaloJetObj &a, const l1CaloJetObj &b) { + return a.jetClusterET > b.jetClusterET; + }); + + /************************************************************************** + * Progress to adding L1EGs built from ECAL TPs 9x9 grid. + * Recall, for 9x9 trigger towers gives diameter 0.78 + ******************************************************************************/ + + // Cluster together the L1EGs around existing HCAL Jet + // Cluster within dEta/dPhi 0.4 which is very close to 0.39 = 9x9/2 + //std::cout << " - Input L1EGs: " << crystalClustersVect.size() << std::endl; + for (auto &caloJetObj : l1CaloJetObjs) { + params["seed_pt"] = caloJetObj.seedTowerET; + params["seed_eta"] = caloJetObj.seedTower.eta(); + params["seed_phi"] = caloJetObj.seedTower.phi(); + params["seed_iEta"] = caloJetObj.seed_iEta; + params["seed_iPhi"] = caloJetObj.seed_iPhi; + params["seed_energy"] = caloJetObj.seedTower.energy(); + + params["hcal_pt"] = caloJetObj.hcalJetClusterET; + params["hcal_seed"] = caloJetObj.hcal_seed; + //params["hcal_3x3"] = caloJetObj.hcal_3x3; + params["hcal_3x5"] = caloJetObj.hcal_3x5; + //params["hcal_5x5"] = caloJetObj.hcal_5x5; + //params["hcal_5x7"] = caloJetObj.hcal_5x7; + params["hcal_7x7"] = caloJetObj.hcal_7x7; + //params["hcal_2x3"] = std::max( caloJetObj.hcal_2x3, std::max( caloJetObj.hcal_2x3_1, caloJetObj.hcal_2x3_2 )); + //params["hcal_2x2"] = std::max( caloJetObj.hcal_2x2, std::max( caloJetObj.hcal_2x2_1, std::max( caloJetObj.hcal_2x2_2, std::max( caloJetObj.hcal_2x2_3, caloJetObj.hcal_2x2_4 )))); + params["hcal_nHits"] = caloJetObj.hcal_nHits; + + params["ecal_pt"] = caloJetObj.ecalJetClusterET; + params["ecal_seed"] = caloJetObj.ecal_seed; + //params["ecal_3x3"] = caloJetObj.ecal_3x3; + params["ecal_3x5"] = caloJetObj.ecal_3x5; + //params["ecal_5x5"] = caloJetObj.ecal_5x5; + //params["ecal_5x7"] = caloJetObj.ecal_5x7; + params["ecal_7x7"] = caloJetObj.ecal_7x7; + //params["ecal_2x3"] = std::max( caloJetObj.ecal_2x3, std::max( caloJetObj.ecal_2x3_1, caloJetObj.ecal_2x3_2 )); + //params["ecal_2x2"] = std::max( caloJetObj.ecal_2x2, std::max( caloJetObj.ecal_2x2_1, std::max( caloJetObj.ecal_2x2_2, std::max( caloJetObj.ecal_2x2_3, caloJetObj.ecal_2x2_4 )))); + params["ecal_nHits"] = caloJetObj.ecal_nHits; + + params["l1eg_pt"] = caloJetObj.l1egJetClusterET; + params["l1eg_seed"] = caloJetObj.l1eg_seed; + //params["l1eg_3x3"] = caloJetObj.l1eg_3x3; + params["l1eg_3x5"] = caloJetObj.l1eg_3x5; + //params["l1eg_5x5"] = caloJetObj.l1eg_5x5; + //params["l1eg_5x7"] = caloJetObj.l1eg_5x7; + params["l1eg_7x7"] = caloJetObj.l1eg_7x7; + //params["l1eg_2x3"] = std::max( caloJetObj.l1eg_2x3, std::max( caloJetObj.l1eg_2x3_1, caloJetObj.l1eg_2x3_2 )); + //params["l1eg_2x2"] = std::max( caloJetObj.l1eg_2x2, std::max( caloJetObj.l1eg_2x2_1, std::max( caloJetObj.l1eg_2x2_2, std::max( caloJetObj.l1eg_2x2_3, caloJetObj.l1eg_2x2_4 )))); + params["l1eg_nHits"] = caloJetObj.l1eg_nHits; + params["l1eg_nL1EGs"] = caloJetObj.l1eg_nL1EGs; + params["l1eg_nL1EGs_standaloneSS"] = caloJetObj.l1eg_nL1EGs_standaloneSS; + params["l1eg_nL1EGs_standaloneIso"] = caloJetObj.l1eg_nL1EGs_standaloneIso; + params["l1eg_nL1EGs_trkMatchSS"] = caloJetObj.l1eg_nL1EGs_trkMatchSS; + params["l1eg_nL1EGs_trkMatchIso"] = caloJetObj.l1eg_nL1EGs_trkMatchIso; + + params["total_et"] = caloJetObj.jetClusterET; + params["total_seed"] = caloJetObj.total_seed; + //params["total_3x3"] = caloJetObj.total_3x3; + params["total_3x5"] = caloJetObj.total_3x5; + //params["total_5x5"] = caloJetObj.total_5x5; + //params["total_5x7"] = caloJetObj.total_5x7; + params["total_7x7"] = caloJetObj.total_7x7; + //params["total_2x3"] = std::max( caloJetObj.total_2x3, std::max( caloJetObj.total_2x3_1, caloJetObj.total_2x3_2 )); + //params["total_2x2"] = std::max( caloJetObj.total_2x2, std::max( caloJetObj.total_2x2_1, std::max( caloJetObj.total_2x2_2, std::max( caloJetObj.total_2x2_3, caloJetObj.total_2x2_4 )))); + params["total_nHits"] = caloJetObj.total_nHits; + //params["total_nTowers"] = total_nTowers; + + //// return -9 for energy and dR values for ecalJet as defaults + float hovere = -9; + if (caloJetObj.ecalJetClusterET > 0.0) { + hovere = caloJetObj.hcalJetClusterET / (caloJetObj.ecalJetClusterET + caloJetObj.l1egJetClusterET); + } + + params["jet_pt"] = caloJetObj.jetClusterET; + params["jet_eta"] = caloJetObj.jetCluster.eta(); + params["jet_phi"] = caloJetObj.jetCluster.phi(); + params["jet_mass"] = caloJetObj.jetCluster.mass(); + params["jet_energy"] = caloJetObj.jetCluster.energy(); + + // Calibrations + params["hcal_calibration"] = + get_hcal_calibration(params["jet_pt"], params["ecal_pt"], params["l1eg_pt"], params["jet_eta"]); + params["hcal_pt_calibration"] = params["hcal_pt"] * params["hcal_calibration"]; + params["jet_pt_calibration"] = params["hcal_pt_calibration"] + params["ecal_pt"] + params["l1eg_pt"]; + + // Tau Vars + // The tau pT calibration is applied as a SF to the total raw pT + // in contrast to the jet calibration above + params["tau_pt_calibration_value"] = get_tau_pt_calibration(params["total_3x5"], + params["ecal_3x5"], + params["l1eg_3x5"], + caloJetObj.n_l1eg_HoverE_LessThreshold, + params["jet_eta"]); + params["tau_pt"] = params["total_3x5"] * params["tau_pt_calibration_value"]; + params["n_l1eg_HoverE_LessThreshold"] = caloJetObj.n_l1eg_HoverE_LessThreshold; + // Currently, applying the tau_pt calibration to the isolation region as well... + // One could switch to using the calibrated jet_pt instead for the iso region... + // This should be revisited - FIXME? + params["tau_total_iso_et"] = params["jet_pt"] * params["tau_pt_calibration_value"]; + params["tau_iso_et"] = (params["jet_pt"] * params["tau_pt_calibration_value"]) - params["tau_pt"]; + params["loose_iso_tau_wp"] = float(loose_iso_tau_wp(params["tau_pt"], params["tau_iso_et"], params["jet_eta"])); + + float calibratedPt = -1; + float ECalIsolation = -1; // Need to loop over 7x7 crystals of unclustered energy + float totalPtPUcorr = -1; + l1tp2::CaloJet caloJet(caloJetObj.jetCluster, calibratedPt, hovere, ECalIsolation, totalPtPUcorr); + caloJet.setExperimentalParams(params); + caloJet.setAssociated_l1EGs(caloJetObj.associated_l1EGs_); + + // Only store jets passing ET threshold + if (params["jet_pt_calibration"] >= EtMinForCollection) { + L1CaloJetsNoCuts->push_back(caloJet); + //L1CaloJetsWithCuts->push_back( caloJet ); + reco::Candidate::PolarLorentzVector jet_p4 = reco::Candidate::PolarLorentzVector( + params["jet_pt_calibration"], caloJet.p4().eta(), caloJet.p4().phi(), caloJet.p4().M()); + L1CaloJetCollectionBXV->push_back(0, l1t::Jet(jet_p4)); + + if (debug) + LogDebug("L1CaloJetProducer") << " Made a Jet, eta " << caloJetObj.jetCluster.eta() << " phi " + << caloJetObj.jetCluster.phi() << " pt " << caloJetObj.jetClusterET + << " calibrated pt " << params["jet_pt_calibration"] << "\n"; + } + + // Only store taus passing ET threshold + if (params["tau_pt"] >= EtMinForTauCollection) { + short int tau_ieta = caloJetObj.seed_iEta; + short int tau_iphi = caloJetObj.seed_iPhi; + short int raw_et = params["total_3x5"]; + short int iso_et = params["tau_iso_et"]; + bool hasEM = false; + if (params["l1eg_3x5"] > 0. || params["ecal_3x5"] > 0.) { + hasEM = true; + } + int tau_qual = int(params["loose_iso_tau_wp"]); + + reco::Candidate::PolarLorentzVector tau_p4 = reco::Candidate::PolarLorentzVector( + params["tau_pt"], caloJet.p4().eta(), caloJet.p4().phi(), caloJet.p4().M()); + l1t::Tau l1Tau = l1t::Tau(tau_p4, params["tau_pt"], caloJet.p4().eta(), caloJet.p4().phi(), tau_qual, iso_et); + l1Tau.setTowerIEta(tau_ieta); + l1Tau.setTowerIPhi(tau_iphi); + l1Tau.setRawEt(raw_et); + l1Tau.setIsoEt(iso_et); + l1Tau.setHasEM(hasEM); + l1Tau.setIsMerged(false); + L1CaloTauCollectionBXV->push_back(0, l1Tau); + + if (debug) { + LogDebug("L1CaloJetProducer") << " Made a Jet, eta " << l1Tau.eta() << " phi " << l1Tau.phi() << " pt " + << l1Tau.rawEt() << " calibrated pt " << l1Tau.pt() << "\n"; + } + } + + } // end jetClusters loop + + iEvent.put(std::move(L1CaloJetsNoCuts), "L1CaloJetsNoCuts"); + //iEvent.put(std::move(L1CaloJetsWithCuts), "L1CaloJetsWithCuts" ); + //iEvent.put(std::move(L1CaloClusterCollectionWithCuts), "L1CaloClusterCollectionWithCuts" ); + iEvent.put(std::move(L1CaloJetCollectionBXV), "L1CaloJetCollectionBXV"); + iEvent.put(std::move(L1CaloTauCollectionBXV), "L1CaloTauCollectionBXV"); +} + +int L1CaloJetProducer::tower_diPhi(int &iPhi_1, int &iPhi_2) const { + // 360 Crystals in full, 72 towers, half way is 36 + int PI = 36; + int result = iPhi_1 - iPhi_2; + while (result > PI) + result -= 2 * PI; + while (result <= -PI) + result += 2 * PI; + return result; +} + +// Added b/c of the iEta jump from +1 to -1 across the barrel mid point +int L1CaloJetProducer::tower_diEta(int &iEta_1, int &iEta_2) const { + // On same side of barrel + if (iEta_1 * iEta_2 > 0) + return iEta_1 - iEta_2; + else + return iEta_1 - iEta_2 - 1; +} + +float L1CaloJetProducer::get_deltaR(reco::Candidate::PolarLorentzVector &p4_1, + reco::Candidate::PolarLorentzVector &p4_2) const { + // Check that pt is > 0 for both or else reco::deltaR returns bogus values + if (p4_1.pt() > 0 && p4_2.pt() > 0) { + return reco::deltaR(p4_1, p4_2); + } else + return -1; +} + +// Apply calibrations to HCAL energy based on EM Fraction, Jet Eta, Jet pT +float L1CaloJetProducer::get_hcal_calibration(float &jet_pt, + float &ecal_pt, + float &ecal_L1EG_jet_pt, + float &jet_eta) const { + float em_frac = (ecal_L1EG_jet_pt + ecal_pt) / jet_pt; + float abs_eta = fabs(jet_eta); + float tmp_jet_pt = jet_pt; + if (tmp_jet_pt > 499) + tmp_jet_pt = 499; + + // Different indices sizes in different calo regions. + // Barrel... + size_t em_index = 0; + size_t eta_index = 0; + size_t pt_index = 0; + float calib = 1.0; + if (abs_eta <= 1.5) { + // Start loop checking 2nd value + for (unsigned int i = 1; i < emFractionBinsBarrel.size(); i++) { + if (em_frac <= emFractionBinsBarrel.at(i)) + break; + em_index++; + } + + // Start loop checking 2nd value + for (unsigned int i = 1; i < absEtaBinsBarrel.size(); i++) { + if (abs_eta <= absEtaBinsBarrel.at(i)) + break; + eta_index++; + } + + // Start loop checking 2nd value + for (unsigned int i = 1; i < jetPtBins.size(); i++) { + if (tmp_jet_pt <= jetPtBins.at(i)) + break; + pt_index++; + } + calib = calibrationsBarrel[eta_index][em_index][pt_index]; + } // end Barrel + else if (abs_eta <= 3.0) // HGCal + { + // Start loop checking 2nd value + for (unsigned int i = 1; i < emFractionBinsHGCal.size(); i++) { + if (em_frac <= emFractionBinsHGCal.at(i)) + break; + em_index++; + } + + // Start loop checking 2nd value + for (unsigned int i = 1; i < absEtaBinsHGCal.size(); i++) { + if (abs_eta <= absEtaBinsHGCal.at(i)) + break; + eta_index++; + } + + // Start loop checking 2nd value + for (unsigned int i = 1; i < jetPtBins.size(); i++) { + if (tmp_jet_pt <= jetPtBins.at(i)) + break; + pt_index++; + } + calib = calibrationsHGCal[eta_index][em_index][pt_index]; + } // end HGCal + else // HF + { + // Start loop checking 2nd value + for (unsigned int i = 1; i < emFractionBinsHF.size(); i++) { + if (em_frac <= emFractionBinsHF.at(i)) + break; + em_index++; + } + + // Start loop checking 2nd value + for (unsigned int i = 1; i < absEtaBinsHF.size(); i++) { + if (abs_eta <= absEtaBinsHF.at(i)) + break; + eta_index++; + } + + // Start loop checking 2nd value + for (unsigned int i = 1; i < jetPtBins.size(); i++) { + if (tmp_jet_pt <= jetPtBins.at(i)) + break; + pt_index++; + } + calib = calibrationsHF[eta_index][em_index][pt_index]; + } // end HF + + return calib; +} + +// Apply calibrations to tau pT based on L1EG info, EM Fraction, Tau Eta, Tau pT +float L1CaloJetProducer::get_tau_pt_calibration( + float &tau_pt, float &ecal_pt, float &l1EG_pt, float &n_L1EGs, float &tau_eta) const { + float em_frac = (l1EG_pt + ecal_pt) / tau_pt; + float abs_eta = fabs(tau_eta); + float tmp_tau_pt = tau_pt; + if (tmp_tau_pt > 199) + tmp_tau_pt = 199; + + // Different indices sizes in different calo regions. + // Barrel... + size_t em_index = 0; + size_t eta_index = 0; + size_t n_L1EG_index = 0; + size_t pt_index = 0; + float calib = 1.0; + // HERE + if (abs_eta <= 1.5) { + // Start loop checking 1st value + for (unsigned int i = 0; i < tauL1egValuesBarrel.size(); i++) { + if (n_L1EGs == tauL1egValuesBarrel.at(i)) + break; + if (tauL1egValuesBarrel.at(i) == tauL1egValuesBarrel.back()) + break; // to preven incrementing on last one + n_L1EG_index++; + } + + // Find key value pair matching n L1EGs + for (auto &l1eg_info : tauL1egInfoMapBarrel) { + if (l1eg_info.first != double(n_L1EG_index)) + continue; + // Start loop checking 2nd value + for (unsigned int i = 1; i < l1eg_info.second.size(); i++) { + if (em_frac <= l1eg_info.second.at(i)) + break; + em_index++; + } + } + + // Start loop checking 2nd value + for (unsigned int i = 1; i < tauAbsEtaBinsBarrel.size(); i++) { + if (abs_eta <= tauAbsEtaBinsBarrel.at(i)) + break; + eta_index++; + } + + // Start loop checking 2nd value + for (unsigned int i = 1; i < tauPtBins.size(); i++) { + if (tmp_tau_pt <= tauPtBins.at(i)) + break; + pt_index++; + } + calib = tauPtCalibrationsBarrel[eta_index][n_L1EG_index][em_index][pt_index]; + } // end Barrel + else if (abs_eta <= 3.0) // HGCal + { + // Start loop checking 1st value + for (unsigned int i = 0; i < tauL1egValuesHGCal.size(); i++) { + if (n_L1EGs == tauL1egValuesHGCal.at(i)) + break; + if (tauL1egValuesHGCal.at(i) == tauL1egValuesHGCal.back()) + break; // to preven incrementing on last one + n_L1EG_index++; + } + + // Find key value pair matching n L1EGs + for (auto &l1eg_info : tauL1egInfoMapHGCal) { + if (l1eg_info.first != double(n_L1EG_index)) + continue; + // Start loop checking 2nd value + for (unsigned int i = 1; i < l1eg_info.second.size(); i++) { + if (em_frac <= l1eg_info.second.at(i)) + break; + em_index++; + } + } + + // Start loop checking 2nd value + for (unsigned int i = 1; i < tauAbsEtaBinsHGCal.size(); i++) { + if (abs_eta <= tauAbsEtaBinsHGCal.at(i)) + break; + eta_index++; + } + + // Start loop checking 2nd value + for (unsigned int i = 1; i < tauPtBins.size(); i++) { + if (tmp_tau_pt <= tauPtBins.at(i)) + break; + pt_index++; + } + calib = tauPtCalibrationsHGCal[eta_index][n_L1EG_index][em_index][pt_index]; + } // end HGCal + else + return calib; + + return calib; +} + +// Loose IsoTau WP +int L1CaloJetProducer::loose_iso_tau_wp(float &tau_pt, float &tau_iso_et, float &tau_eta) const { + // Fully relaxed above 100 GeV pT + if (tau_pt > 100) { + return 1; + } + // Split by barrel and HGCal + // with Barrel first + if (fabs(tau_eta) < 1.5) { + if (isoTauBarrel.Eval(tau_pt) >= (tau_iso_et / tau_pt)) { + return 1; + } else { + return 0; + } + } + // HGCal + if (fabs(tau_eta) < 3.0) { + if (isoTauHGCal.Eval(tau_pt) >= (tau_iso_et / tau_pt)) { + return 1; + } else { + return 0; + } + } + // Beyond HGCal + return 0; +} + +DEFINE_FWK_MODULE(L1CaloJetProducer); diff --git a/L1Trigger/L1CaloTrigger/plugins/L1EGammaCrystalsEmulatorProducer.cc b/L1Trigger/L1CaloTrigger/plugins/L1EGammaCrystalsEmulatorProducer.cc index fc5693620355a..8d8ba03cadec4 100644 --- a/L1Trigger/L1CaloTrigger/plugins/L1EGammaCrystalsEmulatorProducer.cc +++ b/L1Trigger/L1CaloTrigger/plugins/L1EGammaCrystalsEmulatorProducer.cc @@ -240,14 +240,14 @@ class L1EGCrystalClusterEmulatorProducer : public edm::stream::EDProducer<> { edm::EDGetTokenT ecalTPEBToken_; edm::EDGetTokenT > hcalTPToken_; - edm::ESHandle decoder_; + edm::ESGetToken decoderTag_; l1tp2::ParametricCalibration calib_; - edm::ESHandle caloGeometry_; + edm::ESGetToken caloGeometryTag_; const CaloSubdetectorGeometry* ebGeometry; const CaloSubdetectorGeometry* hbGeometry; - edm::ESHandle hbTopology; + edm::ESGetToken hbTopologyTag_; const HcalTopology* hcTopology_; struct mycluster { @@ -341,10 +341,13 @@ L1EGCrystalClusterEmulatorProducer::L1EGCrystalClusterEmulatorProducer(const edm : ecalTPEBToken_(consumes(iConfig.getParameter("ecalTPEB"))), hcalTPToken_( consumes >(iConfig.getParameter("hcalTP"))), - calib_(iConfig.getParameter("calib")) { + decoderTag_(esConsumes(edm::ESInputTag("", ""))), + calib_(iConfig.getParameter("calib")), + caloGeometryTag_(esConsumes(edm::ESInputTag("", ""))), + hbTopologyTag_(esConsumes(edm::ESInputTag("", ""))) { produces(); produces >(); - produces(); + produces("L1CaloTowerCollection"); } L1EGCrystalClusterEmulatorProducer::~L1EGCrystalClusterEmulatorProducer() {} @@ -355,13 +358,14 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: edm::Handle pcalohits; iEvent.getByToken(ecalTPEBToken_, pcalohits); - iSetup.get().get(caloGeometry_); - ebGeometry = caloGeometry_->getSubdetectorGeometry(DetId::Ecal, EcalBarrel); - hbGeometry = caloGeometry_->getSubdetectorGeometry(DetId::Hcal, HcalBarrel); - iSetup.get().get(hbTopology); - hcTopology_ = hbTopology.product(); + const auto& caloGeometry = iSetup.getData(caloGeometryTag_); + ebGeometry = caloGeometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel); + hbGeometry = caloGeometry.getSubdetectorGeometry(DetId::Hcal, HcalBarrel); + const auto& hbTopology = iSetup.getData(hbTopologyTag_); + hcTopology_ = &hbTopology; HcalTrigTowerGeometry theTrigTowerGeometry(hcTopology_); - iSetup.get().get(decoder_); + + const auto& decoder = iSetup.getData(decoderTag_); //**************************************************************** //******************* Get all the hits *************************** @@ -395,7 +399,7 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: edm::Handle > hbhecoll; iEvent.getByToken(hcalTPToken_, hbhecoll); for (const auto& hit : *hbhecoll.product()) { - float et = decoder_->hcaletValue(hit.id(), hit.t0()); + float et = decoder.hcaletValue(hit.id(), hit.t0()); if (et <= 0) continue; if (!(hcTopology_->validHT(hit.id()))) { @@ -1173,7 +1177,7 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: iEvent.put(std::move(L1EGXtalClusters)); iEvent.put(std::move(L1EGammas)); - iEvent.put(std::move(L1CaloTowers)); + iEvent.put(std::move(L1CaloTowers), "L1CaloTowerCollection"); } bool L1EGCrystalClusterEmulatorProducer::passes_iso(float pt, float iso) { diff --git a/L1Trigger/L1CaloTrigger/plugins/L1TowerCalibrator.cc b/L1Trigger/L1CaloTrigger/plugins/L1TowerCalibrator.cc new file mode 100644 index 0000000000000..9054acca9739e --- /dev/null +++ b/L1Trigger/L1CaloTrigger/plugins/L1TowerCalibrator.cc @@ -0,0 +1,579 @@ +// -*- C++ -*- +// +// Package: L1CaloTrigger +// Class: L1TowerCalibrator +// +/**\class L1TowerCalibrator L1TowerCalibrator.cc + +Description: +Take the calibrated unclustered ECAL energy and total HCAL +energy associated with the L1CaloTower collection output +from L1EGammaCrystalsEmulatorProducer: l1CaloTowerCollection, "L1CaloTowerCollection" + +as well as HGCal Tower level inputs: +BXVector "hgcalTriggerPrimitiveDigiProducer" "tower" "HLT" + +and HCAL HF inputs from: +edm::SortedCollection > "simHcalTriggerPrimitiveDigis" "" "HLT" + + +Implement PU-based calibrations which scale down the ET +in the towers based on mapping nTowers with ECAL(HCAL) ET <= defined PU threshold. +This value has been shown to be similar between TTbar, QCD, and minBias samples. +This allows a prediction of nvtx. Which can be mapped to the total minBias +energy in an eta slice of the detector. Subtract the total estimated minBias energy +per eta slice divided by nTowers in that eta slice from each trigger tower in +that eta slice. + +This is all ECAL / HCAL specific or EM vs. Hadronic for HGCal. + +Implementation: +[Notes on implementation] +*/ +// +// Original Author: Tyler Ruggles +// Created: Thu Nov 15 2018 +// $Id$ +// +// + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" + +#include + +#include "DataFormats/L1TCalorimeterPhase2/interface/CaloTower.h" +#include "DataFormats/L1THGCal/interface/HGCalTower.h" +#include "DataFormats/HcalDigi/interface/HcalDigiCollections.h" +#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h" + +#include "CalibFormats/CaloTPG/interface/CaloTPGTranscoder.h" +#include "CalibFormats/CaloTPG/interface/CaloTPGRecord.h" + +#include "L1Trigger/L1TCalorimeter/interface/CaloTools.h" + +#include "TFile.h" +#include "TF1.h" + +class L1TowerCalibrator : public edm::EDProducer { +public: + explicit L1TowerCalibrator(const edm::ParameterSet&); + +private: + void produce(edm::Event&, const edm::EventSetup&) override; + + double HcalTpEtMin; + double EcalTpEtMin; + double HGCalHadTpEtMin; + double HGCalEmTpEtMin; + double HFTpEtMin; + double puThreshold; + double puThresholdL1eg; + double puThresholdHcalMin; + double puThresholdHcalMax; + double puThresholdEcalMin; + double puThresholdEcalMax; + double puThresholdHGCalEMMin; + double puThresholdHGCalEMMax; + double puThresholdHGCalHadMin; + double puThresholdHGCalHadMax; + double puThresholdHFMin; + double puThresholdHFMax; + double barrelSF; + double hgcalSF; + double hfSF; + bool debug; + bool skipCalibrations; + + edm::EDGetTokenT l1TowerToken_; + edm::Handle l1CaloTowerHandle; + + edm::EDGetTokenT hgcalTowersToken_; + edm::Handle hgcalTowersHandle; + l1t::HGCalTowerBxCollection hgcalTowers; + + edm::EDGetTokenT hcalToken_; + edm::Handle hcalTowerHandle; + edm::ESGetToken decoderTag_; + + // nHits to nvtx functions + std::vector nHits_to_nvtx_params; + std::map nHits_to_nvtx_funcs; + + // nvtx to PU subtraction functions + std::vector nvtx_to_PU_sub_params; + std::map ecal_nvtx_to_PU_sub_funcs; + std::map hcal_nvtx_to_PU_sub_funcs; + std::map hgcalEM_nvtx_to_PU_sub_funcs; + std::map hgcalHad_nvtx_to_PU_sub_funcs; + std::map hf_nvtx_to_PU_sub_funcs; + std::map > all_nvtx_to_PU_sub_funcs; +}; + +L1TowerCalibrator::L1TowerCalibrator(const edm::ParameterSet& iConfig) + : HcalTpEtMin(iConfig.getParameter("HcalTpEtMin")), + EcalTpEtMin(iConfig.getParameter("EcalTpEtMin")), + HGCalHadTpEtMin(iConfig.getParameter("HGCalHadTpEtMin")), + HGCalEmTpEtMin(iConfig.getParameter("HGCalEmTpEtMin")), + HFTpEtMin(iConfig.getParameter("HFTpEtMin")), + puThreshold(iConfig.getParameter("puThreshold")), + puThresholdL1eg(iConfig.getParameter("puThresholdL1eg")), + puThresholdHcalMin(iConfig.getParameter("puThresholdHcalMin")), + puThresholdHcalMax(iConfig.getParameter("puThresholdHcalMax")), + puThresholdEcalMin(iConfig.getParameter("puThresholdEcalMin")), + puThresholdEcalMax(iConfig.getParameter("puThresholdEcalMax")), + puThresholdHGCalEMMin(iConfig.getParameter("puThresholdHGCalEMMin")), + puThresholdHGCalEMMax(iConfig.getParameter("puThresholdHGCalEMMax")), + puThresholdHGCalHadMin(iConfig.getParameter("puThresholdHGCalHadMin")), + puThresholdHGCalHadMax(iConfig.getParameter("puThresholdHGCalHadMax")), + puThresholdHFMin(iConfig.getParameter("puThresholdHFMin")), + puThresholdHFMax(iConfig.getParameter("puThresholdHFMax")), + barrelSF(iConfig.getParameter("barrelSF")), + hgcalSF(iConfig.getParameter("hgcalSF")), + hfSF(iConfig.getParameter("hfSF")), + debug(iConfig.getParameter("debug")), + skipCalibrations(iConfig.getParameter("skipCalibrations")), + l1TowerToken_(consumes(iConfig.getParameter("l1CaloTowers"))), + hgcalTowersToken_( + consumes(iConfig.getParameter("L1HgcalTowersInputTag"))), + hcalToken_(consumes(iConfig.getParameter("hcalDigis"))), + decoderTag_(esConsumes(edm::ESInputTag("", ""))), + nHits_to_nvtx_params(iConfig.getParameter >("nHits_to_nvtx_params")), + nvtx_to_PU_sub_params(iConfig.getParameter >("nvtx_to_PU_sub_params")) { + // Initialize the nHits --> nvtx functions + for (uint i = 0; i < nHits_to_nvtx_params.size(); i++) { + edm::ParameterSet* pset = &nHits_to_nvtx_params.at(i); + std::string calo = pset->getParameter("fit"); + nHits_to_nvtx_funcs[calo.c_str()] = TF1(calo.c_str(), "[0] + [1] * x"); + nHits_to_nvtx_funcs[calo.c_str()].SetParameter(0, pset->getParameter >("params").at(0)); + nHits_to_nvtx_funcs[calo.c_str()].SetParameter(1, pset->getParameter >("params").at(1)); + + if (debug) { + printf( + "nHits_to_nvtx_params[%i]\n \ + fit: %s \n \ + p1: %f \n \ + p2 %f \n", + i, + calo.c_str(), + nHits_to_nvtx_funcs[calo.c_str()].GetParameter(0), + nHits_to_nvtx_funcs[calo.c_str()].GetParameter(1)); + } + } + + // Initialize the nvtx --> PU subtraction functions + all_nvtx_to_PU_sub_funcs["ecal"] = ecal_nvtx_to_PU_sub_funcs; + all_nvtx_to_PU_sub_funcs["hcal"] = hcal_nvtx_to_PU_sub_funcs; + all_nvtx_to_PU_sub_funcs["hgcalEM"] = hgcalEM_nvtx_to_PU_sub_funcs; + all_nvtx_to_PU_sub_funcs["hgcalHad"] = hgcalHad_nvtx_to_PU_sub_funcs; + all_nvtx_to_PU_sub_funcs["hf"] = hf_nvtx_to_PU_sub_funcs; + + for (uint i = 0; i < nvtx_to_PU_sub_params.size(); i++) { + edm::ParameterSet* pset = &nvtx_to_PU_sub_params.at(i); + std::string calo = pset->getParameter("calo"); + std::string iEta = pset->getParameter("iEta"); + double p1 = pset->getParameter >("params").at(0); + double p2 = pset->getParameter >("params").at(1); + + all_nvtx_to_PU_sub_funcs[calo.c_str()][iEta.c_str()] = TF1(calo.c_str(), "[0] + [1] * x"); + all_nvtx_to_PU_sub_funcs[calo.c_str()][iEta.c_str()].SetParameter(0, p1); + all_nvtx_to_PU_sub_funcs[calo.c_str()][iEta.c_str()].SetParameter(1, p2); + + if (debug) { + printf( + "nvtx_to_PU_sub_params[%i]\n \ + sub detector: %s \n \ + iEta: %s \n \ + p1: %f \n \ + p2 %f \n", + i, + calo.c_str(), + iEta.c_str(), + p1, + p2); + } + } + + // Our two outputs, calibrated towers and estimated nvtx for fun + produces("L1CaloTowerCalibratedCollection"); + produces("EstimatedNvtx"); +} + +void L1TowerCalibrator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + // Estimated number of vertices used for calibration estimattion + std::unique_ptr EstimatedNvtx(new double); + // Calibrated output collection + std::unique_ptr L1CaloTowerCalibratedCollection(new l1tp2::CaloTowerCollection); + + // Load the ECAL+HCAL tower sums coming from L1EGammaCrystalsEmulatorProducer.cc + iEvent.getByToken(l1TowerToken_, l1CaloTowerHandle); + + // HGCal info + iEvent.getByToken(hgcalTowersToken_, hgcalTowersHandle); + hgcalTowers = (*hgcalTowersHandle.product()); + + // HF Tower info + iEvent.getByToken(hcalToken_, hcalTowerHandle); + + // Barrel ECAL (unclustered) and HCAL + for (auto& hit : *l1CaloTowerHandle.product()) { + l1tp2::CaloTower l1Hit; + l1Hit.setEcalTowerEt(hit.ecalTowerEt()); + l1Hit.setHcalTowerEt(hit.hcalTowerEt()); + l1Hit.setL1egTowerEt(hit.l1egTowerEt()); + // Add min ET thresholds for tower ET + if (l1Hit.ecalTowerEt() < EcalTpEtMin) + l1Hit.setEcalTowerEt(0.0); + if (l1Hit.hcalTowerEt() < HcalTpEtMin) + l1Hit.setHcalTowerEt(0.0); + l1Hit.setTowerIEta(hit.towerIEta()); + l1Hit.setTowerIPhi(hit.towerIPhi()); + l1Hit.setTowerEta(hit.towerEta()); + l1Hit.setTowerPhi(hit.towerPhi()); + l1Hit.setIsBarrel(hit.isBarrel()); + l1Hit.setNL1eg(hit.nL1eg()); + l1Hit.setL1egTrkSS(hit.l1egTrkSS()); + l1Hit.setL1egTrkIso(hit.l1egTrkIso()); + l1Hit.setL1egStandaloneSS(hit.l1egStandaloneSS()); + l1Hit.setL1egStandaloneIso(hit.l1egStandaloneIso()); + + // FIXME There is an error in the L1EGammaCrystalsEmulatorProducer.cc which is + // returning towers with minimal ECAL energy, and no HCAL energy with these + // iEta/iPhi coordinates and eta = -88.653152 and phi = -99.000000. + // Skip these for the time being until the upstream code has been debugged + if ((int)l1Hit.towerIEta() == -1016 && (int)l1Hit.towerIPhi() == -962) + continue; + + (*L1CaloTowerCalibratedCollection).push_back(l1Hit); + if (debug) + printf("Barrel tower iEta %i iPhi %i eta %f phi %f ecal_et %f hcal_et_sum %f\n", + (int)l1Hit.towerIEta(), + (int)l1Hit.towerIPhi(), + l1Hit.towerEta(), + l1Hit.towerPhi(), + l1Hit.ecalTowerEt(), + l1Hit.hcalTowerEt()); + } + + // Loop over HGCalTowers and create L1CaloTowers for them and add to collection + // This iterator is taken from the PF P2 group + // https://github.com/p2l1pfp/cmssw/blob/170808db68038d53794bc65fdc962f8fc337a24d/L1Trigger/Phase2L1ParticleFlow/plugins/L1TPFCaloProducer.cc#L278-L289 + for (auto it = hgcalTowers.begin(0), ed = hgcalTowers.end(0); it != ed; ++it) { + // skip lowest ET towers + if (it->etEm() < HGCalEmTpEtMin && it->etHad() < HGCalHadTpEtMin) + continue; + + l1tp2::CaloTower l1Hit; + // Set energies normally, but need to zero if below threshold + if (it->etEm() < HGCalEmTpEtMin) + l1Hit.setEcalTowerEt(0.); + else + l1Hit.setEcalTowerEt(it->etEm()); + + if (it->etHad() < HGCalHadTpEtMin) + l1Hit.setHcalTowerEt(0.); + else + l1Hit.setHcalTowerEt(it->etHad()); + + l1Hit.setTowerEta(it->eta()); + l1Hit.setTowerPhi(it->phi()); + l1Hit.setTowerIEta(-98); // -98 mean HGCal + l1Hit.setTowerIPhi(-98); + l1Hit.setIsBarrel(false); + (*L1CaloTowerCalibratedCollection).push_back(l1Hit); + if (debug) + printf("HGCal tower iEta %i iPhi %i eta %f phi %f ecal_et %f hcal_et_sum %f\n", + (int)l1Hit.towerIEta(), + (int)l1Hit.towerIPhi(), + l1Hit.towerEta(), + l1Hit.towerPhi(), + l1Hit.ecalTowerEt(), + l1Hit.hcalTowerEt()); + } + + // Loop over Hcal HF tower inputs and create L1CaloTowers and add to + // L1CaloTowerCalibratedCollection collection + const auto& decoder = iSetup.getData(decoderTag_); + for (const auto& hit : *hcalTowerHandle.product()) { + HcalTrigTowerDetId id = hit.id(); + double et = decoder.hcaletValue(hit.id(), hit.t0()); + if (et < HFTpEtMin) + continue; + // Only doing HF so skip outside range + if (abs(id.ieta()) < l1t::CaloTools::kHFBegin) + continue; + if (abs(id.ieta()) > l1t::CaloTools::kHFEnd) + continue; + + l1tp2::CaloTower l1Hit; + l1Hit.setEcalTowerEt(0.); + l1Hit.setHcalTowerEt(et); + l1Hit.setTowerEta(l1t::CaloTools::towerEta(id.ieta())); + l1Hit.setTowerPhi(l1t::CaloTools::towerPhi(id.ieta(), id.iphi())); + l1Hit.setTowerIEta(id.ieta()); + l1Hit.setTowerIPhi(id.iphi()); + l1Hit.setIsBarrel(false); + (*L1CaloTowerCalibratedCollection).push_back(l1Hit); + + if (debug) + printf("HCAL HF tower iEta %i iPhi %i eta %f phi %f ecal_et %f hcal_et_sum %f\n", + (int)l1Hit.towerIEta(), + (int)l1Hit.towerIPhi(), + l1Hit.towerEta(), + l1Hit.towerPhi(), + l1Hit.ecalTowerEt(), + l1Hit.hcalTowerEt()); + } + + // N Tower totals + // For mapping to estimated nvtx in event + int i_ecal_hits_leq_threshold = 0; + int i_hgcalEM_hits_leq_threshold = 0; + int i_hcal_hits_leq_threshold = 0; + int i_hgcalHad_hits_leq_threshold = 0; + int i_hf_hits_leq_threshold = 0; + + // Loop over the collection containing all hits + // and calculate the number of hits falling into the + // "less than or equal" nTowers variable which maps to + // estimated number of vertices + for (auto& l1CaloTower : (*L1CaloTowerCalibratedCollection)) { + // Barrel ECAL + if (l1CaloTower.ecalTowerEt() > 0. && l1CaloTower.towerIEta() != -98) { + if (l1CaloTower.ecalTowerEt() <= puThresholdEcalMax && l1CaloTower.ecalTowerEt() >= puThresholdEcalMin) { + i_ecal_hits_leq_threshold++; + } + } + + // HGCal EM + if (l1CaloTower.ecalTowerEt() > 0. && l1CaloTower.towerIEta() == -98) { + if (l1CaloTower.ecalTowerEt() <= puThresholdHGCalEMMax && l1CaloTower.ecalTowerEt() >= puThresholdHGCalEMMin) { + i_hgcalEM_hits_leq_threshold++; + } + } + + // Barrel HCAL + if (l1CaloTower.hcalTowerEt() > 0. && l1CaloTower.towerIEta() != -98 && + abs(l1CaloTower.towerEta()) < 2.0) // abs(eta) < 2 just keeps us out of HF + { + if (l1CaloTower.hcalTowerEt() <= puThresholdHcalMax && l1CaloTower.hcalTowerEt() >= puThresholdHcalMin) { + i_hcal_hits_leq_threshold++; + } + } + + // HGCal Had + if (l1CaloTower.hcalTowerEt() > 0. && l1CaloTower.towerIEta() == -98) { + if (l1CaloTower.hcalTowerEt() <= puThresholdHGCalHadMax && l1CaloTower.hcalTowerEt() >= puThresholdHGCalHadMin) { + i_hgcalHad_hits_leq_threshold++; + } + } + + // HF + if (l1CaloTower.hcalTowerEt() > 0. && l1CaloTower.towerIEta() != -98 && + abs(l1CaloTower.towerEta()) > 2.0) // abs(eta) > 2 keeps us out of barrel HF + { + if (l1CaloTower.hcalTowerEt() <= puThresholdHFMax && l1CaloTower.hcalTowerEt() >= puThresholdHFMin) { + i_hf_hits_leq_threshold++; + } + } + } + + // For each subdetector, map to the estimated number of vertices + double ecal_nvtx = nHits_to_nvtx_funcs["ecal"].Eval(i_ecal_hits_leq_threshold); + double hcal_nvtx = nHits_to_nvtx_funcs["hcal"].Eval(i_hcal_hits_leq_threshold); + double hgcalEM_nvtx = nHits_to_nvtx_funcs["hgcalEM"].Eval(i_hgcalEM_hits_leq_threshold); + double hgcalHad_nvtx = nHits_to_nvtx_funcs["hgcalHad"].Eval(i_hgcalHad_hits_leq_threshold); + double hf_nvtx = nHits_to_nvtx_funcs["hf"].Eval(i_hf_hits_leq_threshold); + // Make sure all values are >= 0 + if (ecal_nvtx < 0) + ecal_nvtx = 0; + if (hcal_nvtx < 0) + hcal_nvtx = 0; + if (hgcalEM_nvtx < 0) + hgcalEM_nvtx = 0; + if (hgcalHad_nvtx < 0) + hgcalHad_nvtx = 0; + if (hf_nvtx < 0) + hf_nvtx = 0; + // Best estimate is avg of all except HF. + // This is b/c HF currently has such poor prediction power, it only degrades the avg result + // NEW, with 10_3_X, hgcal and HF has the best results based on the values I took... + // skip ECAL and HCAL + //*EstimatedNvtx = ( ecal_nvtx + hcal_nvtx + hgcalEM_nvtx + hgcalHad_nvtx + hf_nvtx ) / 3.; + *EstimatedNvtx = (hgcalEM_nvtx + hgcalHad_nvtx + hf_nvtx) / 3.; + + if (debug) { + double lumi = iEvent.eventAuxiliary().luminosityBlock(); + double event = iEvent.eventAuxiliary().event(); + + printf( + "L1TowerCalibrater: lumi %.0f evt %.0f nTowers for subdetecters \ + \nECAL: %i --> nvtx = %.1f \ + \nHGCal EM: %i --> nvtx = %.1f \ + \nHCAL: %i --> nvtx = %.1f \ + \nHGCal Had: %i --> nvtx = %.1f \ + \nHCAL HF: %i --> nvtx = %.1f \ + \nEstimated Nvtx = %.1f\n", + lumi, + event, + i_ecal_hits_leq_threshold, + ecal_nvtx, + i_hgcalEM_hits_leq_threshold, + hgcalEM_nvtx, + i_hcal_hits_leq_threshold, + hcal_nvtx, + i_hgcalHad_hits_leq_threshold, + hgcalHad_nvtx, + i_hf_hits_leq_threshold, + hf_nvtx, + *EstimatedNvtx); + } + + // Use estimated number of vertices to subtract off PU contributions + // to each and every hit. In cases where the energy would go negative, + // limit this to zero. + if (!skipCalibrations) // skipCalibrations simply passes the towers through + { + for (auto& l1CaloTower : (*L1CaloTowerCalibratedCollection)) { + // Barrel ECAL eta slices + if (l1CaloTower.ecalTowerEt() > 0. && l1CaloTower.towerIEta() != -98) { + if (abs(l1CaloTower.towerIEta()) <= 3) { + l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() - + all_nvtx_to_PU_sub_funcs["ecal"]["er1to3"].Eval(*EstimatedNvtx) * barrelSF); + } + if (abs(l1CaloTower.towerIEta()) <= 6 && abs(l1CaloTower.towerIEta()) >= 4) { + l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() - + all_nvtx_to_PU_sub_funcs["ecal"]["er4to6"].Eval(*EstimatedNvtx) * barrelSF); + } + if (abs(l1CaloTower.towerIEta()) <= 9 && abs(l1CaloTower.towerIEta()) >= 7) { + l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() - + all_nvtx_to_PU_sub_funcs["ecal"]["er7to9"].Eval(*EstimatedNvtx) * barrelSF); + } + if (abs(l1CaloTower.towerIEta()) <= 12 && abs(l1CaloTower.towerIEta()) >= 10) { + l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() - + all_nvtx_to_PU_sub_funcs["ecal"]["er10to12"].Eval(*EstimatedNvtx) * barrelSF); + } + if (abs(l1CaloTower.towerIEta()) <= 15 && abs(l1CaloTower.towerIEta()) >= 13) { + l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() - + all_nvtx_to_PU_sub_funcs["ecal"]["er13to15"].Eval(*EstimatedNvtx) * barrelSF); + } + if (abs(l1CaloTower.towerIEta()) <= 18 && abs(l1CaloTower.towerIEta()) >= 16) { + l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() - + all_nvtx_to_PU_sub_funcs["ecal"]["er16to18"].Eval(*EstimatedNvtx) * barrelSF); + } + } + + // HGCal EM eta slices + if (l1CaloTower.ecalTowerEt() > 0. && l1CaloTower.towerIEta() == -98) { + if (abs(l1CaloTower.towerEta()) <= 1.8) { + l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() - + all_nvtx_to_PU_sub_funcs["hgcalEM"]["er1p4to1p8"].Eval(*EstimatedNvtx) * hgcalSF); + } + if (abs(l1CaloTower.towerEta()) <= 2.1 && abs(l1CaloTower.towerEta()) > 1.8) { + l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() - + all_nvtx_to_PU_sub_funcs["hgcalEM"]["er1p8to2p1"].Eval(*EstimatedNvtx) * hgcalSF); + } + if (abs(l1CaloTower.towerEta()) <= 2.4 && abs(l1CaloTower.towerEta()) > 2.1) { + l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() - + all_nvtx_to_PU_sub_funcs["hgcalEM"]["er2p1to2p4"].Eval(*EstimatedNvtx) * hgcalSF); + } + if (abs(l1CaloTower.towerEta()) <= 2.7 && abs(l1CaloTower.towerEta()) > 2.4) { + l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() - + all_nvtx_to_PU_sub_funcs["hgcalEM"]["er2p4to2p7"].Eval(*EstimatedNvtx) * hgcalSF); + } + if (abs(l1CaloTower.towerEta()) <= 3.1 && abs(l1CaloTower.towerEta()) > 2.7) { + l1CaloTower.setEcalTowerEt(l1CaloTower.ecalTowerEt() - + all_nvtx_to_PU_sub_funcs["hgcalEM"]["er2p7to3p1"].Eval(*EstimatedNvtx) * hgcalSF); + } + } + + // Barrel HCAL eta slices + if (l1CaloTower.hcalTowerEt() > 0. && l1CaloTower.towerIEta() != -98 && + abs(l1CaloTower.towerEta()) < 2.0) // abs(eta) < 2 just keeps us out of HF + { + if (abs(l1CaloTower.towerIEta()) <= 3) { + l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() - + all_nvtx_to_PU_sub_funcs["hcal"]["er1to3"].Eval(*EstimatedNvtx) * barrelSF); + } + if (abs(l1CaloTower.towerIEta()) <= 6 && abs(l1CaloTower.towerIEta()) >= 4) { + l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() - + all_nvtx_to_PU_sub_funcs["hcal"]["er4to6"].Eval(*EstimatedNvtx) * barrelSF); + } + if (abs(l1CaloTower.towerIEta()) <= 9 && abs(l1CaloTower.towerIEta()) >= 7) { + l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() - + all_nvtx_to_PU_sub_funcs["hcal"]["er7to9"].Eval(*EstimatedNvtx) * barrelSF); + } + if (abs(l1CaloTower.towerIEta()) <= 12 && abs(l1CaloTower.towerIEta()) >= 10) { + l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() - + all_nvtx_to_PU_sub_funcs["hcal"]["er10to12"].Eval(*EstimatedNvtx) * barrelSF); + } + if (abs(l1CaloTower.towerIEta()) <= 15 && abs(l1CaloTower.towerIEta()) >= 13) { + l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() - + all_nvtx_to_PU_sub_funcs["hcal"]["er13to15"].Eval(*EstimatedNvtx) * barrelSF); + } + if (abs(l1CaloTower.towerIEta()) <= 18 && abs(l1CaloTower.towerIEta()) >= 16) { + l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() - + all_nvtx_to_PU_sub_funcs["hcal"]["er16to18"].Eval(*EstimatedNvtx) * barrelSF); + } + } + + // HGCal Had eta slices + if (l1CaloTower.hcalTowerEt() > 0. && l1CaloTower.towerIEta() == -98) { + if (abs(l1CaloTower.towerEta()) <= 1.8) { + l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() - + all_nvtx_to_PU_sub_funcs["hgcalHad"]["er1p4to1p8"].Eval(*EstimatedNvtx) * hgcalSF); + } + if (abs(l1CaloTower.towerEta()) <= 2.1 && abs(l1CaloTower.towerEta()) > 1.8) { + l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() - + all_nvtx_to_PU_sub_funcs["hgcalHad"]["er1p8to2p1"].Eval(*EstimatedNvtx) * hgcalSF); + } + if (abs(l1CaloTower.towerEta()) <= 2.4 && abs(l1CaloTower.towerEta()) > 2.1) { + l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() - + all_nvtx_to_PU_sub_funcs["hgcalHad"]["er2p1to2p4"].Eval(*EstimatedNvtx) * hgcalSF); + } + if (abs(l1CaloTower.towerEta()) <= 2.7 && abs(l1CaloTower.towerEta()) > 2.4) { + l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() - + all_nvtx_to_PU_sub_funcs["hgcalHad"]["er2p4to2p7"].Eval(*EstimatedNvtx) * hgcalSF); + } + if (abs(l1CaloTower.towerEta()) <= 3.1 && abs(l1CaloTower.towerEta()) > 2.7) { + l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() - + all_nvtx_to_PU_sub_funcs["hgcalHad"]["er2p7to3p1"].Eval(*EstimatedNvtx) * hgcalSF); + } + } + + // HF eta slices + if (l1CaloTower.hcalTowerEt() > 0. && l1CaloTower.towerIEta() != -98 && + abs(l1CaloTower.towerEta()) > 2.0) // abs(eta) > 2 keeps us out of barrel HF + { + if (abs(l1CaloTower.towerIEta()) <= 33 && abs(l1CaloTower.towerIEta()) >= 29) { + l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() - + all_nvtx_to_PU_sub_funcs["hf"]["er29to33"].Eval(*EstimatedNvtx) * hfSF); + } + if (abs(l1CaloTower.towerIEta()) <= 37 && abs(l1CaloTower.towerIEta()) >= 34) { + l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() - + all_nvtx_to_PU_sub_funcs["hf"]["er34to37"].Eval(*EstimatedNvtx) * hfSF); + } + if (abs(l1CaloTower.towerIEta()) <= 41 && abs(l1CaloTower.towerIEta()) >= 38) { + l1CaloTower.setHcalTowerEt(l1CaloTower.hcalTowerEt() - + all_nvtx_to_PU_sub_funcs["hf"]["er38to41"].Eval(*EstimatedNvtx) * hfSF); + } + } + + // Make sure none are negative + if (l1CaloTower.ecalTowerEt() < 0.) + l1CaloTower.setEcalTowerEt(0.); + if (l1CaloTower.hcalTowerEt() < 0.) + l1CaloTower.setHcalTowerEt(0.); + } + } + + iEvent.put(std::move(EstimatedNvtx), "EstimatedNvtx"); + iEvent.put(std::move(L1CaloTowerCalibratedCollection), "L1CaloTowerCalibratedCollection"); +} + +DEFINE_FWK_MODULE(L1TowerCalibrator); diff --git a/L1Trigger/L1CaloTrigger/python/L1CaloJetHTTProducer_cfi.py b/L1Trigger/L1CaloTrigger/python/L1CaloJetHTTProducer_cfi.py new file mode 100644 index 0000000000000..11650c362daf3 --- /dev/null +++ b/L1Trigger/L1CaloTrigger/python/L1CaloJetHTTProducer_cfi.py @@ -0,0 +1,10 @@ +import FWCore.ParameterSet.Config as cms + +L1CaloJetHTTProducer = cms.EDProducer("L1CaloJetHTTProducer", + EtaMax = cms.double(2.4), + PtMin = cms.double(30.0), + BXVCaloJetsInputTag = cms.InputTag("L1CaloJetProducer","L1CaloJetCollectionBXV"), + genJets = cms.InputTag("ak4GenJetsNoNu", "", "HLT"), + debug = cms.bool(False), + use_gen_jets = cms.bool(False), +) diff --git a/L1Trigger/L1CaloTrigger/python/L1CaloJetProducer_cfi.py b/L1Trigger/L1CaloTrigger/python/L1CaloJetProducer_cfi.py new file mode 100644 index 0000000000000..88be72570b512 --- /dev/null +++ b/L1Trigger/L1CaloTrigger/python/L1CaloJetProducer_cfi.py @@ -0,0 +1,164 @@ +import FWCore.ParameterSet.Config as cms + +L1CaloJetProducer = cms.EDProducer("L1CaloJetProducer", + debug = cms.bool(False), + HcalTpEtMin = cms.double(0.5), + EcalTpEtMin = cms.double(0.5), + HGCalHadTpEtMin = cms.double(0.5), + HGCalEmTpEtMin = cms.double(0.5), + HFTpEtMin = cms.double(0.5), + EtMinForSeedHit = cms.double(2.5), + EtMinForCollection = cms.double(10), + EtMinForTauCollection = cms.double(10), + l1CaloTowers = cms.InputTag("L1TowerCalibrationProducer","L1CaloTowerCalibratedCollection"), + L1CrystalClustersInputTag = cms.InputTag("L1EGammaClusterEmuProducer", "L1EGXtalClusterEmulator"), + #L1HgcalTowersInputTag = cms.InputTag("hgcalTriggerPrimitiveDigiProducer","tower"), + #hcalDigis = cms.InputTag("simHcalTriggerPrimitiveDigis"), + + # Calibrations derived 18 April 2019 on 10_3_X MTD QCD sample + jetPtBins = cms.vdouble([ 0.0,5.0,7.5,10.0,12.5,15.0,17.5,20.0,22.5,25.0,27.5,30.0,35.0,40.0,45.0,50.0,55.0,60.0,65.0,70.0,75.0,80.0,85.0,90.0,95.0,100.0,110.0,120.0,130.0,140.0,150.0,160.0,170.0,180.0,190.0,200.0,225.0,250.0,275.0,300.0,325.0,400.0,500.0]), + emFractionBinsBarrel = cms.vdouble([ 0.00,0.31,0.40,0.47,0.53,0.58,0.63,0.69,0.76,0.84,1.05]), + absEtaBinsBarrel = cms.vdouble([ 0.00,0.30,0.60,1.00,1.50]), + jetCalibrationsBarrel = cms.vdouble([ + 1.640, 1.626, 1.617, 1.608, 1.599, 1.591, 1.583, 1.574, 1.566, 1.559, 1.551, 1.540, 1.525, 1.512, 1.498, 1.486, 1.474, 1.462, 1.451, 1.441, 1.431, 1.421, 1.412, 1.404, 1.396, 1.385, 1.372, 1.360, 1.350, 1.341, 1.334, 1.329, 1.325, 1.322, 1.320, 1.320, 1.325, 1.336, 1.352, 1.373, 1.428, 1.553, + 1.855, 1.839, 1.829, 1.819, 1.809, 1.800, 1.790, 1.781, 1.771, 1.762, 1.753, 1.740, 1.723, 1.706, 1.690, 1.674, 1.659, 1.645, 1.631, 1.617, 1.604, 1.591, 1.579, 1.568, 1.556, 1.540, 1.520, 1.501, 1.484, 1.469, 1.454, 1.441, 1.430, 1.419, 1.410, 1.396, 1.382, 1.374, 1.370, 1.372, 1.386, 1.443, + 1.993, 1.978, 1.967, 1.957, 1.947, 1.937, 1.927, 1.918, 1.908, 1.898, 1.889, 1.875, 1.857, 1.839, 1.822, 1.805, 1.789, 1.773, 1.757, 1.742, 1.727, 1.713, 1.699, 1.685, 1.672, 1.653, 1.629, 1.606, 1.584, 1.564, 1.545, 1.527, 1.511, 1.495, 1.481, 1.459, 1.434, 1.414, 1.400, 1.392, 1.391, 1.431, + 2.155, 2.135, 2.123, 2.110, 2.097, 2.085, 2.073, 2.061, 2.049, 2.037, 2.026, 2.009, 1.987, 1.965, 1.944, 1.924, 1.905, 1.885, 1.867, 1.849, 1.832, 1.815, 1.798, 1.783, 1.767, 1.745, 1.717, 1.691, 1.666, 1.643, 1.622, 1.602, 1.583, 1.566, 1.550, 1.525, 1.495, 1.470, 1.451, 1.437, 1.422, 1.424, + 2.305, 2.284, 2.271, 2.257, 2.244, 2.231, 2.218, 2.205, 2.193, 2.180, 2.168, 2.150, 2.126, 2.103, 2.081, 2.059, 2.037, 2.017, 1.996, 1.976, 1.957, 1.938, 1.920, 1.902, 1.884, 1.859, 1.827, 1.797, 1.768, 1.741, 1.715, 1.691, 1.668, 1.647, 1.627, 1.595, 1.556, 1.523, 1.497, 1.477, 1.451, 1.446, + 2.555, 2.528, 2.510, 2.492, 2.474, 2.457, 2.440, 2.423, 2.406, 2.390, 2.373, 2.349, 2.318, 2.288, 2.258, 2.229, 2.201, 2.174, 2.148, 2.122, 2.097, 2.072, 2.049, 2.026, 2.004, 1.971, 1.930, 1.892, 1.856, 1.822, 1.790, 1.761, 1.733, 1.707, 1.683, 1.645, 1.600, 1.563, 1.535, 1.515, 1.493, 1.505, + 2.850, 2.815, 2.792, 2.770, 2.747, 2.725, 2.703, 2.682, 2.661, 2.640, 2.619, 2.589, 2.549, 2.511, 2.473, 2.437, 2.401, 2.367, 2.334, 2.301, 2.270, 2.239, 2.209, 2.181, 2.153, 2.112, 2.062, 2.014, 1.969, 1.927, 1.888, 1.852, 1.818, 1.787, 1.758, 1.712, 1.658, 1.616, 1.585, 1.563, 1.544, 1.576, + 3.386, 3.325, 3.286, 3.248, 3.210, 3.174, 3.138, 3.103, 3.069, 3.036, 3.004, 2.956, 2.896, 2.838, 2.783, 2.730, 2.680, 2.632, 2.586, 2.543, 2.501, 2.461, 2.423, 2.387, 2.352, 2.303, 2.242, 2.187, 2.136, 2.090, 2.047, 2.008, 1.973, 1.940, 1.910, 1.863, 1.806, 1.760, 1.721, 1.688, 1.636, 1.572, + 4.450, 4.344, 4.276, 4.210, 4.146, 4.083, 4.022, 3.962, 3.904, 3.847, 3.792, 3.712, 3.610, 3.513, 3.420, 3.333, 3.250, 3.171, 3.096, 3.025, 2.958, 2.894, 2.833, 2.776, 2.721, 2.645, 2.552, 2.468, 2.393, 2.325, 2.265, 2.210, 2.162, 2.118, 2.079, 2.021, 1.956, 1.908, 1.873, 1.849, 1.822, 1.818, + 6.912, 6.764, 6.667, 6.571, 6.477, 6.384, 6.293, 6.203, 6.114, 6.027, 5.941, 5.814, 5.650, 5.491, 5.336, 5.187, 5.042, 4.903, 4.767, 4.637, 4.511, 4.389, 4.272, 4.159, 4.050, 3.894, 3.700, 3.522, 3.359, 3.210, 3.076, 2.954, 2.846, 2.750, 2.667, 2.549, 2.438, 2.390, 2.400, 2.463, 2.731, 3.577, + 1.536, 1.527, 1.521, 1.515, 1.509, 1.503, 1.497, 1.491, 1.486, 1.480, 1.475, 1.467, 1.457, 1.447, 1.437, 1.428, 1.419, 1.411, 1.403, 1.395, 1.388, 1.380, 1.374, 1.367, 1.361, 1.353, 1.342, 1.333, 1.325, 1.319, 1.313, 1.309, 1.305, 1.303, 1.302, 1.302, 1.307, 1.317, 1.333, 1.354, 1.409, 1.544, + 1.740, 1.730, 1.723, 1.716, 1.709, 1.702, 1.696, 1.689, 1.683, 1.676, 1.670, 1.661, 1.648, 1.636, 1.624, 1.613, 1.602, 1.590, 1.580, 1.569, 1.559, 1.549, 1.539, 1.530, 1.521, 1.507, 1.490, 1.474, 1.459, 1.445, 1.432, 1.420, 1.409, 1.399, 1.389, 1.375, 1.360, 1.350, 1.345, 1.346, 1.363, 1.439, + 1.869, 1.858, 1.851, 1.843, 1.836, 1.829, 1.822, 1.815, 1.808, 1.801, 1.795, 1.785, 1.771, 1.758, 1.746, 1.733, 1.721, 1.709, 1.697, 1.686, 1.674, 1.663, 1.652, 1.642, 1.631, 1.616, 1.597, 1.578, 1.560, 1.543, 1.527, 1.512, 1.498, 1.485, 1.472, 1.452, 1.427, 1.408, 1.393, 1.383, 1.377, 1.407, + 2.010, 1.997, 1.988, 1.980, 1.971, 1.963, 1.954, 1.946, 1.938, 1.930, 1.922, 1.910, 1.894, 1.879, 1.864, 1.849, 1.834, 1.820, 1.806, 1.792, 1.779, 1.766, 1.753, 1.740, 1.728, 1.709, 1.686, 1.664, 1.642, 1.622, 1.603, 1.584, 1.567, 1.550, 1.534, 1.509, 1.478, 1.451, 1.430, 1.414, 1.396, 1.407, + 2.171, 2.155, 2.144, 2.134, 2.123, 2.113, 2.103, 2.093, 2.083, 2.073, 2.063, 2.049, 2.030, 2.012, 1.994, 1.976, 1.959, 1.942, 1.925, 1.909, 1.893, 1.877, 1.862, 1.847, 1.833, 1.812, 1.785, 1.759, 1.734, 1.710, 1.688, 1.667, 1.646, 1.627, 1.609, 1.579, 1.542, 1.509, 1.482, 1.459, 1.426, 1.401, + 2.366, 2.345, 2.332, 2.319, 2.306, 2.293, 2.280, 2.268, 2.255, 2.243, 2.231, 2.213, 2.189, 2.166, 2.144, 2.122, 2.100, 2.079, 2.058, 2.038, 2.018, 1.999, 1.980, 1.962, 1.944, 1.918, 1.884, 1.853, 1.822, 1.794, 1.766, 1.741, 1.716, 1.693, 1.671, 1.636, 1.591, 1.554, 1.522, 1.497, 1.461, 1.441, + 2.680, 2.647, 2.626, 2.605, 2.584, 2.564, 2.544, 2.525, 2.506, 2.487, 2.469, 2.442, 2.407, 2.373, 2.340, 2.309, 2.278, 2.249, 2.221, 2.193, 2.167, 2.142, 2.117, 2.093, 2.070, 2.038, 1.996, 1.958, 1.922, 1.888, 1.857, 1.828, 1.801, 1.775, 1.751, 1.713, 1.666, 1.626, 1.592, 1.563, 1.516, 1.459, + 3.093, 3.051, 3.023, 2.995, 2.969, 2.942, 2.916, 2.890, 2.865, 2.840, 2.816, 2.780, 2.733, 2.688, 2.644, 2.602, 2.561, 2.522, 2.484, 2.447, 2.411, 2.376, 2.343, 2.311, 2.280, 2.235, 2.178, 2.126, 2.077, 2.031, 1.988, 1.949, 1.912, 1.878, 1.846, 1.796, 1.735, 1.685, 1.646, 1.614, 1.571, 1.543, + 4.080, 3.996, 3.941, 3.888, 3.836, 3.785, 3.735, 3.686, 3.638, 3.591, 3.546, 3.479, 3.393, 3.312, 3.233, 3.159, 3.087, 3.019, 2.953, 2.891, 2.831, 2.774, 2.720, 2.668, 2.618, 2.548, 2.462, 2.383, 2.312, 2.247, 2.188, 2.134, 2.086, 2.042, 2.003, 1.943, 1.875, 1.825, 1.788, 1.763, 1.738, 1.747, + 6.659, 6.494, 6.387, 6.282, 6.179, 6.078, 5.978, 5.881, 5.785, 5.691, 5.599, 5.463, 5.289, 5.122, 4.961, 4.806, 4.657, 4.514, 4.378, 4.246, 4.121, 4.001, 3.886, 3.776, 3.671, 3.522, 3.341, 3.177, 3.030, 2.899, 2.783, 2.682, 2.594, 2.519, 2.456, 2.373, 2.311, 2.307, 2.355, 2.447, 2.743, 3.536, + 1.519, 1.512, 1.507, 1.502, 1.498, 1.493, 1.488, 1.484, 1.479, 1.475, 1.471, 1.464, 1.456, 1.448, 1.440, 1.433, 1.425, 1.418, 1.412, 1.405, 1.399, 1.393, 1.387, 1.381, 1.376, 1.369, 1.359, 1.351, 1.344, 1.337, 1.332, 1.327, 1.323, 1.320, 1.318, 1.317, 1.319, 1.326, 1.339, 1.356, 1.404, 1.530, + 1.719, 1.710, 1.705, 1.699, 1.693, 1.688, 1.682, 1.677, 1.671, 1.666, 1.661, 1.653, 1.642, 1.632, 1.622, 1.613, 1.603, 1.594, 1.584, 1.575, 1.567, 1.558, 1.550, 1.541, 1.533, 1.522, 1.507, 1.492, 1.479, 1.466, 1.454, 1.443, 1.432, 1.422, 1.413, 1.399, 1.382, 1.369, 1.361, 1.356, 1.359, 1.401, + 1.860, 1.851, 1.844, 1.838, 1.832, 1.826, 1.820, 1.814, 1.808, 1.802, 1.796, 1.787, 1.775, 1.764, 1.753, 1.742, 1.731, 1.720, 1.710, 1.699, 1.689, 1.679, 1.669, 1.660, 1.650, 1.636, 1.618, 1.601, 1.584, 1.568, 1.553, 1.538, 1.524, 1.511, 1.498, 1.478, 1.451, 1.428, 1.410, 1.394, 1.374, 1.373, + 1.993, 1.982, 1.974, 1.967, 1.960, 1.953, 1.946, 1.938, 1.931, 1.924, 1.917, 1.907, 1.894, 1.880, 1.867, 1.854, 1.841, 1.828, 1.816, 1.804, 1.792, 1.780, 1.768, 1.757, 1.746, 1.729, 1.708, 1.687, 1.667, 1.648, 1.629, 1.612, 1.595, 1.578, 1.563, 1.537, 1.504, 1.475, 1.451, 1.430, 1.402, 1.388, + 2.146, 2.131, 2.122, 2.113, 2.104, 2.095, 2.086, 2.077, 2.068, 2.059, 2.051, 2.038, 2.021, 2.004, 1.988, 1.972, 1.957, 1.941, 1.926, 1.911, 1.897, 1.883, 1.869, 1.855, 1.842, 1.822, 1.797, 1.773, 1.749, 1.727, 1.706, 1.685, 1.665, 1.647, 1.629, 1.599, 1.561, 1.527, 1.497, 1.472, 1.431, 1.389, + 2.336, 2.319, 2.307, 2.296, 2.285, 2.273, 2.262, 2.251, 2.240, 2.230, 2.219, 2.203, 2.182, 2.162, 2.142, 2.122, 2.103, 2.084, 2.065, 2.047, 2.029, 2.012, 1.994, 1.978, 1.961, 1.937, 1.906, 1.876, 1.848, 1.820, 1.794, 1.769, 1.745, 1.722, 1.700, 1.665, 1.618, 1.578, 1.543, 1.512, 1.465, 1.420, + 2.594, 2.571, 2.556, 2.541, 2.526, 2.511, 2.497, 2.482, 2.468, 2.453, 2.439, 2.419, 2.391, 2.364, 2.338, 2.312, 2.287, 2.262, 2.238, 2.214, 2.191, 2.168, 2.146, 2.124, 2.103, 2.072, 2.032, 1.993, 1.957, 1.922, 1.889, 1.857, 1.827, 1.799, 1.772, 1.728, 1.672, 1.625, 1.585, 1.552, 1.507, 1.482, + 3.016, 2.982, 2.959, 2.937, 2.915, 2.893, 2.872, 2.851, 2.830, 2.809, 2.788, 2.758, 2.718, 2.679, 2.641, 2.604, 2.568, 2.533, 2.499, 2.465, 2.432, 2.400, 2.369, 2.339, 2.309, 2.266, 2.211, 2.160, 2.111, 2.064, 2.020, 1.979, 1.940, 1.904, 1.870, 1.815, 1.749, 1.694, 1.651, 1.618, 1.580, 1.590, + 4.054, 3.966, 3.909, 3.854, 3.800, 3.747, 3.696, 3.646, 3.597, 3.549, 3.503, 3.435, 3.348, 3.266, 3.187, 3.112, 3.041, 2.974, 2.909, 2.848, 2.790, 2.734, 2.681, 2.631, 2.583, 2.516, 2.434, 2.360, 2.292, 2.232, 2.178, 2.128, 2.084, 2.045, 2.009, 1.955, 1.895, 1.850, 1.817, 1.794, 1.769, 1.766, + 6.752, 6.570, 6.452, 6.337, 6.224, 6.113, 6.004, 5.898, 5.795, 5.693, 5.594, 5.449, 5.263, 5.085, 4.915, 4.752, 4.597, 4.449, 4.308, 4.173, 4.045, 3.923, 3.808, 3.697, 3.593, 3.446, 3.269, 3.110, 2.969, 2.846, 2.738, 2.644, 2.564, 2.497, 2.442, 2.372, 2.323, 2.327, 2.374, 2.459, 2.718, 3.371, + 1.615, 1.605, 1.599, 1.593, 1.587, 1.581, 1.575, 1.570, 1.564, 1.558, 1.553, 1.545, 1.534, 1.524, 1.514, 1.504, 1.494, 1.485, 1.476, 1.468, 1.460, 1.452, 1.444, 1.437, 1.429, 1.419, 1.407, 1.395, 1.385, 1.375, 1.367, 1.360, 1.353, 1.348, 1.344, 1.339, 1.336, 1.340, 1.349, 1.364, 1.410, 1.537, + 1.867, 1.856, 1.849, 1.841, 1.834, 1.827, 1.820, 1.813, 1.806, 1.799, 1.792, 1.782, 1.768, 1.755, 1.742, 1.730, 1.718, 1.705, 1.694, 1.682, 1.671, 1.660, 1.649, 1.638, 1.628, 1.613, 1.594, 1.575, 1.558, 1.542, 1.526, 1.512, 1.498, 1.485, 1.474, 1.455, 1.433, 1.417, 1.406, 1.401, 1.404, 1.458, + 2.037, 2.025, 2.017, 2.009, 2.001, 1.993, 1.985, 1.977, 1.970, 1.962, 1.954, 1.943, 1.928, 1.914, 1.899, 1.885, 1.871, 1.858, 1.844, 1.831, 1.818, 1.805, 1.793, 1.781, 1.769, 1.751, 1.728, 1.707, 1.686, 1.666, 1.647, 1.629, 1.612, 1.596, 1.580, 1.555, 1.524, 1.498, 1.477, 1.462, 1.445, 1.461, + 2.206, 2.192, 2.183, 2.174, 2.165, 2.156, 2.148, 2.139, 2.130, 2.122, 2.113, 2.101, 2.084, 2.068, 2.052, 2.037, 2.021, 2.006, 1.992, 1.977, 1.963, 1.949, 1.935, 1.922, 1.908, 1.889, 1.864, 1.840, 1.817, 1.795, 1.773, 1.753, 1.733, 1.714, 1.696, 1.666, 1.628, 1.593, 1.562, 1.536, 1.492, 1.444, + 2.401, 2.385, 2.375, 2.365, 2.355, 2.346, 2.336, 2.326, 2.317, 2.307, 2.298, 2.284, 2.265, 2.247, 2.229, 2.212, 2.195, 2.178, 2.161, 2.144, 2.128, 2.112, 2.097, 2.081, 2.066, 2.044, 2.015, 1.988, 1.961, 1.935, 1.910, 1.886, 1.862, 1.840, 1.818, 1.782, 1.735, 1.692, 1.653, 1.618, 1.559, 1.487, + 2.637, 2.620, 2.609, 2.599, 2.588, 2.577, 2.567, 2.556, 2.546, 2.536, 2.525, 2.510, 2.490, 2.470, 2.451, 2.432, 2.413, 2.394, 2.376, 2.357, 2.339, 2.322, 2.304, 2.287, 2.270, 2.245, 2.213, 2.182, 2.151, 2.122, 2.093, 2.065, 2.038, 2.012, 1.986, 1.944, 1.887, 1.835, 1.787, 1.742, 1.665, 1.561, + 2.952, 2.935, 2.924, 2.913, 2.902, 2.891, 2.879, 2.868, 2.858, 2.847, 2.836, 2.820, 2.798, 2.777, 2.756, 2.736, 2.715, 2.695, 2.675, 2.655, 2.635, 2.616, 2.597, 2.578, 2.559, 2.531, 2.495, 2.459, 2.424, 2.390, 2.357, 2.325, 2.293, 2.262, 2.233, 2.182, 2.114, 2.050, 1.992, 1.938, 1.844, 1.721, + 3.403, 3.388, 3.378, 3.368, 3.359, 3.349, 3.339, 3.330, 3.320, 3.311, 3.301, 3.287, 3.268, 3.249, 3.231, 3.213, 3.194, 3.176, 3.158, 3.140, 3.123, 3.105, 3.088, 3.070, 3.053, 3.028, 2.994, 2.961, 2.928, 2.896, 2.865, 2.834, 2.803, 2.773, 2.743, 2.692, 2.621, 2.554, 2.488, 2.425, 2.307, 2.120, + 4.386, 4.363, 4.347, 4.332, 4.317, 4.302, 4.287, 4.272, 4.258, 4.243, 4.229, 4.208, 4.180, 4.153, 4.127, 4.101, 4.075, 4.051, 4.026, 4.003, 3.979, 3.956, 3.934, 3.912, 3.891, 3.860, 3.821, 3.783, 3.747, 3.712, 3.680, 3.649, 3.620, 3.592, 3.566, 3.524, 3.471, 3.427, 3.390, 3.361, 3.323, 3.312, + 6.918, 6.864, 6.829, 6.794, 6.759, 6.724, 6.690, 6.656, 6.623, 6.589, 6.556, 6.507, 6.443, 6.380, 6.319, 6.258, 6.199, 6.142, 6.085, 6.030, 5.976, 5.923, 5.872, 5.821, 5.773, 5.702, 5.612, 5.528, 5.449, 5.376, 5.308, 5.245, 5.189, 5.138, 5.092, 5.027, 4.964, 4.939, 4.953, 5.005, 5.233, 6.041, + ]), + emFractionBinsHGCal = cms.vdouble([ 0.00,0.55,0.67,0.77,1.05]), + absEtaBinsHGCal = cms.vdouble([ 1.50,1.90,2.40,3.00]), + jetCalibrationsHGCal = cms.vdouble([ + 1.395, 1.394, 1.394, 1.394, 1.394, 1.394, 1.393, 1.393, 1.393, 1.393, 1.393, 1.393, 1.392, 1.392, 1.392, 1.391, 1.391, 1.391, 1.390, 1.390, 1.390, 1.389, 1.389, 1.389, 1.388, 1.388, 1.387, 1.387, 1.386, 1.385, 1.385, 1.384, 1.383, 1.383, 1.382, 1.381, 1.379, 1.378, 1.376, 1.374, 1.371, 1.365, + 1.575, 1.574, 1.574, 1.574, 1.574, 1.574, 1.573, 1.573, 1.573, 1.573, 1.572, 1.572, 1.572, 1.571, 1.571, 1.570, 1.570, 1.569, 1.569, 1.569, 1.568, 1.568, 1.567, 1.567, 1.566, 1.566, 1.565, 1.564, 1.563, 1.562, 1.561, 1.560, 1.559, 1.559, 1.558, 1.556, 1.554, 1.552, 1.549, 1.547, 1.543, 1.535, + 1.846, 1.845, 1.845, 1.844, 1.844, 1.843, 1.843, 1.842, 1.842, 1.841, 1.841, 1.840, 1.839, 1.839, 1.838, 1.837, 1.836, 1.835, 1.834, 1.833, 1.832, 1.831, 1.830, 1.829, 1.829, 1.827, 1.825, 1.824, 1.822, 1.820, 1.818, 1.816, 1.814, 1.813, 1.811, 1.808, 1.803, 1.799, 1.794, 1.790, 1.780, 1.765, + 3.377, 3.355, 3.341, 3.326, 3.312, 3.298, 3.283, 3.269, 3.255, 3.242, 3.228, 3.207, 3.180, 3.154, 3.128, 3.102, 3.077, 3.051, 3.027, 3.002, 2.979, 2.955, 2.932, 2.909, 2.886, 2.853, 2.811, 2.769, 2.729, 2.691, 2.654, 2.618, 2.584, 2.551, 2.519, 2.467, 2.400, 2.341, 2.291, 2.248, 2.186, 2.149, + 1.085, 1.087, 1.088, 1.089, 1.090, 1.091, 1.092, 1.093, 1.094, 1.095, 1.097, 1.098, 1.100, 1.103, 1.105, 1.107, 1.109, 1.111, 1.113, 1.116, 1.118, 1.120, 1.122, 1.124, 1.126, 1.130, 1.134, 1.138, 1.143, 1.147, 1.151, 1.156, 1.160, 1.164, 1.169, 1.176, 1.187, 1.198, 1.209, 1.220, 1.242, 1.280, + 61.559, 11.087, 4.083, 1.971, 1.335, 1.145, 1.090, 1.074, 1.072, 1.072, 1.074, 1.078, 1.083, 1.088, 1.093, 1.098, 1.103, 1.108, 1.112, 1.117, 1.122, 1.127, 1.132, 1.137, 1.142, 1.149, 1.159, 1.169, 1.179, 1.189, 1.199, 1.209, 1.218, 1.228, 1.238, 1.255, 1.280, 1.305, 1.329, 1.354, 1.403, 1.490, + 89.821, 15.567, 5.416, 2.387, 1.485, 1.218, 1.140, 1.119, 1.115, 1.116, 1.118, 1.123, 1.129, 1.135, 1.141, 1.147, 1.153, 1.159, 1.165, 1.171, 1.178, 1.184, 1.190, 1.196, 1.202, 1.211, 1.223, 1.236, 1.248, 1.260, 1.272, 1.284, 1.297, 1.309, 1.321, 1.342, 1.373, 1.403, 1.434, 1.464, 1.525, 1.632, + 83.760, 19.529, 8.102, 3.961, 2.461, 1.918, 1.721, 1.650, 1.624, 1.615, 1.612, 1.611, 1.611, 1.611, 1.612, 1.612, 1.613, 1.613, 1.614, 1.614, 1.615, 1.616, 1.616, 1.617, 1.617, 1.618, 1.619, 1.620, 1.621, 1.622, 1.623, 1.624, 1.625, 1.626, 1.628, 1.629, 1.632, 1.635, 1.637, 1.640, 1.645, 1.655, + 1.262, 1.261, 1.260, 1.259, 1.259, 1.258, 1.257, 1.256, 1.256, 1.255, 1.254, 1.253, 1.252, 1.250, 1.249, 1.247, 1.246, 1.244, 1.243, 1.242, 1.240, 1.239, 1.237, 1.236, 1.234, 1.232, 1.229, 1.226, 1.223, 1.221, 1.218, 1.215, 1.212, 1.209, 1.206, 1.201, 1.194, 1.186, 1.179, 1.172, 1.157, 1.132, + 1.328, 1.327, 1.326, 1.326, 1.325, 1.325, 1.324, 1.324, 1.323, 1.323, 1.322, 1.322, 1.321, 1.320, 1.319, 1.318, 1.317, 1.316, 1.315, 1.314, 1.313, 1.312, 1.311, 1.310, 1.309, 1.308, 1.306, 1.304, 1.302, 1.300, 1.298, 1.296, 1.294, 1.292, 1.290, 1.287, 1.282, 1.277, 1.272, 1.267, 1.257, 1.240, + 1.463, 1.461, 1.460, 1.459, 1.458, 1.457, 1.456, 1.455, 1.454, 1.453, 1.452, 1.451, 1.449, 1.447, 1.445, 1.443, 1.440, 1.438, 1.436, 1.434, 1.432, 1.430, 1.428, 1.426, 1.424, 1.421, 1.417, 1.413, 1.409, 1.405, 1.401, 1.397, 1.393, 1.389, 1.385, 1.377, 1.367, 1.357, 1.347, 1.337, 1.316, 1.281, + 44.899, 14.731, 7.692, 4.568, 3.180, 2.562, 2.285, 2.159, 2.101, 2.072, 2.057, 2.043, 2.031, 2.021, 2.011, 2.001, 1.991, 1.982, 1.972, 1.962, 1.952, 1.942, 1.932, 1.922, 1.912, 1.897, 1.878, 1.858, 1.838, 1.818, 1.798, 1.779, 1.759, 1.739, 1.719, 1.685, 1.635, 1.586, 1.536, 1.487, 1.387, 1.214, + ]), + emFractionBinsHF = cms.vdouble([ 0.00,1.05]), + absEtaBinsHF = cms.vdouble([ 3.00,3.60,6.00]), + jetCalibrationsHF = cms.vdouble([ + 3.223, 2.140, 1.683, 1.364, 1.142, 0.987, 0.879, 0.804, 0.752, 0.716, 0.691, 0.667, 0.651, 0.644, 0.641, 0.641, 0.641, 0.642, 0.644, 0.645, 0.647, 0.648, 0.650, 0.651, 0.653, 0.655, 0.658, 0.662, 0.665, 0.668, 0.671, 0.674, 0.677, 0.681, 0.684, 0.689, 0.697, 0.705, 0.713, 0.721, 0.736, 0.764, + 2.598, 1.813, 1.462, 1.206, 1.019, 0.883, 0.783, 0.711, 0.658, 0.620, 0.592, 0.564, 0.543, 0.532, 0.526, 0.524, 0.523, 0.524, 0.524, 0.525, 0.526, 0.527, 0.528, 0.529, 0.530, 0.532, 0.534, 0.536, 0.538, 0.541, 0.543, 0.545, 0.547, 0.549, 0.552, 0.556, 0.561, 0.567, 0.572, 0.578, 0.589, 0.608, + ]), + + # Testing tau calibrations derived for March workshop, 20 March 2019 + tauPtBins = cms.vdouble([ 0.0,5.0,7.5,10.0,12.5,15.0,20.0,25.0,30.0,35.0,40.0,45.0,50.0,55.0,60.0,70.0,80.0,100.0,150.0,200.0]), + tauAbsEtaBinsBarrel = cms.vdouble([ 0.00,0.30,0.60,1.00,1.50]), + tauL1egInfoBarrel = cms.VPSet( + cms.PSet( + l1egCount = cms.double( 0.0 ), + l1egEmFractions = cms.vdouble([ 0.0, 0.091, 0.317, 1.05]), + ), + cms.PSet( + l1egCount = cms.double( 1.0 ), + l1egEmFractions = cms.vdouble([ 0.0, 0.634, 0.888, 1.05]), + ), + cms.PSet( + l1egCount = cms.double( 2.0 ), + l1egEmFractions = cms.vdouble([ 0.0, 0.821, 0.957, 1.05]), + ), + ), + tauCalibrationsBarrel = cms.vdouble([ + 1.917, 1.714, 1.605, 1.514, 1.437, 1.344, 1.252, 1.188, 1.142, 1.110, 1.087, 1.071, 1.059, 1.051, 1.044, 1.038, 1.034, 1.032, 1.032, + 1.856, 1.729, 1.655, 1.589, 1.530, 1.453, 1.369, 1.301, 1.248, 1.205, 1.170, 1.143, 1.121, 1.103, 1.084, 1.066, 1.050, 1.037, 1.034, + 2.053, 1.884, 1.788, 1.702, 1.627, 1.530, 1.426, 1.345, 1.282, 1.233, 1.194, 1.164, 1.141, 1.122, 1.102, 1.084, 1.070, 1.059, 1.057, + 1.783, 1.659, 1.588, 1.524, 1.466, 1.391, 1.309, 1.244, 1.191, 1.149, 1.116, 1.089, 1.068, 1.050, 1.031, 1.013, 0.998, 0.986, 0.983, + 1.896, 1.713, 1.613, 1.527, 1.453, 1.362, 1.270, 1.202, 1.152, 1.115, 1.088, 1.069, 1.054, 1.044, 1.033, 1.024, 1.018, 1.015, 1.014, + 1.936, 1.677, 1.544, 1.437, 1.350, 1.249, 1.156, 1.095, 1.055, 1.028, 1.011, 1.000, 0.992, 0.987, 0.983, 0.980, 0.978, 0.978, 0.978, + 1.850, 1.683, 1.590, 1.509, 1.439, 1.351, 1.258, 1.188, 1.135, 1.096, 1.065, 1.043, 1.025, 1.012, 0.999, 0.987, 0.979, 0.973, 0.972, + 2.066, 1.784, 1.638, 1.517, 1.418, 1.301, 1.189, 1.114, 1.063, 1.029, 1.006, 0.990, 0.979, 0.972, 0.965, 0.961, 0.958, 0.957, 0.957, + 1.074, 1.063, 1.056, 1.050, 1.043, 1.034, 1.022, 1.011, 1.001, 0.992, 0.983, 0.974, 0.966, 0.959, 0.949, 0.937, 0.921, 0.895, 0.874, + 2.053, 1.757, 1.609, 1.489, 1.394, 1.285, 1.186, 1.123, 1.083, 1.057, 1.040, 1.029, 1.022, 1.018, 1.014, 1.012, 1.010, 1.010, 1.010, + 1.733, 1.628, 1.567, 1.513, 1.463, 1.399, 1.328, 1.270, 1.224, 1.187, 1.157, 1.133, 1.114, 1.098, 1.081, 1.064, 1.050, 1.038, 1.035, + 2.009, 1.852, 1.762, 1.682, 1.611, 1.518, 1.417, 1.337, 1.274, 1.224, 1.184, 1.153, 1.128, 1.108, 1.086, 1.066, 1.049, 1.036, 1.033, + 1.806, 1.669, 1.591, 1.522, 1.461, 1.382, 1.296, 1.230, 1.177, 1.136, 1.103, 1.078, 1.058, 1.043, 1.025, 1.010, 0.998, 0.988, 0.986, + 1.944, 1.725, 1.609, 1.512, 1.431, 1.334, 1.239, 1.173, 1.126, 1.094, 1.071, 1.056, 1.045, 1.037, 1.030, 1.024, 1.021, 1.019, 1.019, + 1.900, 1.631, 1.496, 1.389, 1.304, 1.208, 1.122, 1.068, 1.034, 1.012, 0.999, 0.990, 0.984, 0.981, 0.978, 0.976, 0.975, 0.975, 0.975, + 1.682, 1.577, 1.516, 1.460, 1.409, 1.342, 1.266, 1.203, 1.152, 1.109, 1.075, 1.046, 1.022, 1.003, 0.980, 0.958, 0.937, 0.918, 0.913, + 1.928, 1.693, 1.569, 1.466, 1.381, 1.280, 1.183, 1.116, 1.070, 1.038, 1.016, 1.001, 0.991, 0.984, 0.977, 0.973, 0.970, 0.969, 0.969, + 1.067, 1.058, 1.052, 1.047, 1.041, 1.033, 1.023, 1.014, 1.005, 0.997, 0.989, 0.981, 0.974, 0.968, 0.959, 0.948, 0.933, 0.908, 0.886, + 1.883, 1.663, 1.548, 1.453, 1.375, 1.282, 1.193, 1.133, 1.091, 1.063, 1.044, 1.031, 1.022, 1.016, 1.010, 1.006, 1.004, 1.003, 1.003, + 1.688, 1.595, 1.541, 1.492, 1.448, 1.390, 1.325, 1.271, 1.228, 1.193, 1.164, 1.140, 1.121, 1.105, 1.087, 1.070, 1.055, 1.041, 1.037, + 2.251, 2.020, 1.891, 1.779, 1.682, 1.561, 1.434, 1.339, 1.267, 1.213, 1.173, 1.142, 1.119, 1.102, 1.084, 1.069, 1.058, 1.051, 1.050, + 1.799, 1.654, 1.572, 1.501, 1.438, 1.359, 1.275, 1.211, 1.162, 1.125, 1.096, 1.074, 1.057, 1.044, 1.030, 1.019, 1.009, 1.003, 1.002, + 1.780, 1.630, 1.547, 1.475, 1.413, 1.335, 1.254, 1.193, 1.148, 1.114, 1.089, 1.070, 1.056, 1.045, 1.034, 1.025, 1.018, 1.014, 1.013, + 1.750, 1.549, 1.445, 1.359, 1.288, 1.205, 1.126, 1.072, 1.036, 1.012, 0.995, 0.984, 0.976, 0.971, 0.966, 0.963, 0.961, 0.960, 0.960, + 1.736, 1.610, 1.537, 1.472, 1.414, 1.339, 1.257, 1.192, 1.141, 1.100, 1.067, 1.041, 1.021, 1.005, 0.986, 0.970, 0.956, 0.945, 0.942, + 1.961, 1.712, 1.582, 1.474, 1.386, 1.281, 1.181, 1.114, 1.068, 1.036, 1.015, 1.001, 0.991, 0.984, 0.978, 0.974, 0.971, 0.970, 0.970, + 1.065, 1.056, 1.051, 1.046, 1.040, 1.033, 1.023, 1.014, 1.005, 0.997, 0.989, 0.982, 0.974, 0.968, 0.958, 0.946, 0.930, 0.900, 0.872, + 2.047, 1.810, 1.683, 1.575, 1.485, 1.375, 1.266, 1.188, 1.133, 1.094, 1.066, 1.046, 1.032, 1.022, 1.012, 1.005, 1.000, 0.998, 0.997, + 2.022, 1.861, 1.768, 1.686, 1.613, 1.519, 1.418, 1.338, 1.275, 1.226, 1.187, 1.157, 1.133, 1.114, 1.093, 1.074, 1.059, 1.047, 1.045, + 2.503, 2.236, 2.087, 1.956, 1.843, 1.699, 1.548, 1.434, 1.347, 1.281, 1.231, 1.192, 1.164, 1.142, 1.118, 1.099, 1.084, 1.074, 1.072, + 2.253, 2.007, 1.872, 1.755, 1.654, 1.529, 1.401, 1.305, 1.235, 1.182, 1.143, 1.115, 1.093, 1.078, 1.061, 1.048, 1.039, 1.033, 1.033, + 2.228, 2.010, 1.888, 1.782, 1.690, 1.574, 1.453, 1.361, 1.292, 1.239, 1.200, 1.170, 1.147, 1.130, 1.112, 1.097, 1.085, 1.078, 1.077, + 2.527, 2.136, 1.934, 1.769, 1.633, 1.475, 1.326, 1.226, 1.159, 1.114, 1.084, 1.064, 1.051, 1.042, 1.034, 1.028, 1.025, 1.024, 1.024, + 2.097, 1.908, 1.801, 1.706, 1.622, 1.514, 1.398, 1.308, 1.238, 1.183, 1.140, 1.106, 1.080, 1.060, 1.038, 1.018, 1.002, 0.991, 0.988, + 2.479, 2.106, 1.912, 1.753, 1.623, 1.470, 1.326, 1.229, 1.163, 1.120, 1.090, 1.070, 1.057, 1.048, 1.040, 1.034, 1.031, 1.029, 1.029, + 2.097, 1.906, 1.796, 1.699, 1.612, 1.500, 1.378, 1.281, 1.205, 1.145, 1.098, 1.060, 1.031, 1.007, 0.981, 0.958, 0.939, 0.923, 0.920, + ]), + tauAbsEtaBinsHGCal = cms.vdouble([ 1.50,1.90,2.40,3.00]), + tauL1egInfoHGCal = cms.VPSet( + cms.PSet( + l1egCount = cms.double( 0.0 ), + l1egEmFractions = cms.vdouble([ 0.0, 0.473, 0.72, 0.894, 1.05]), + ), + ), + tauCalibrationsHGCal = cms.vdouble([ + 1.665, 1.556, 1.495, 1.443, 1.398, 1.342, 1.285, 1.242, 1.211, 1.188, 1.170, 1.158, 1.148, 1.141, 1.134, 1.128, 1.124, 1.122, 1.121, + 1.646, 1.525, 1.459, 1.402, 1.353, 1.293, 1.231, 1.186, 1.153, 1.129, 1.112, 1.099, 1.089, 1.082, 1.075, 1.069, 1.065, 1.063, 1.063, + 1.826, 1.630, 1.528, 1.444, 1.376, 1.294, 1.218, 1.166, 1.130, 1.107, 1.091, 1.080, 1.072, 1.067, 1.063, 1.060, 1.058, 1.057, 1.057, + 2.031, 1.762, 1.624, 1.514, 1.424, 1.321, 1.226, 1.164, 1.124, 1.098, 1.081, 1.069, 1.062, 1.057, 1.053, 1.050, 1.049, 1.048, 1.048, + 1.396, 1.328, 1.291, 1.259, 1.232, 1.198, 1.164, 1.139, 1.120, 1.107, 1.097, 1.089, 1.084, 1.080, 1.076, 1.073, 1.071, 1.070, 1.070, + 1.242, 1.214, 1.197, 1.182, 1.168, 1.149, 1.128, 1.111, 1.097, 1.085, 1.075, 1.067, 1.060, 1.055, 1.048, 1.042, 1.036, 1.030, 1.029, + 1.491, 1.368, 1.305, 1.253, 1.211, 1.161, 1.115, 1.084, 1.064, 1.050, 1.041, 1.035, 1.031, 1.028, 1.026, 1.024, 1.023, 1.023, 1.023, + 1.348, 1.289, 1.255, 1.225, 1.199, 1.165, 1.127, 1.098, 1.075, 1.057, 1.043, 1.031, 1.022, 1.016, 1.008, 1.001, 0.995, 0.991, 0.990, + 1.422, 1.373, 1.345, 1.320, 1.297, 1.266, 1.232, 1.205, 1.183, 1.165, 1.150, 1.138, 1.128, 1.120, 1.111, 1.102, 1.095, 1.088, 1.086, + 1.350, 1.322, 1.304, 1.288, 1.272, 1.250, 1.224, 1.201, 1.180, 1.162, 1.146, 1.132, 1.119, 1.108, 1.093, 1.078, 1.061, 1.038, 1.026, + 1.440, 1.373, 1.334, 1.300, 1.270, 1.231, 1.189, 1.155, 1.129, 1.108, 1.092, 1.079, 1.069, 1.061, 1.053, 1.045, 1.038, 1.033, 1.032, + 1.356, 1.317, 1.293, 1.271, 1.250, 1.221, 1.186, 1.156, 1.129, 1.105, 1.085, 1.067, 1.051, 1.037, 1.019, 1.000, 0.980, 0.954, 0.941, + ]), +) + diff --git a/L1Trigger/L1CaloTrigger/python/L1CaloJets_cff.py b/L1Trigger/L1CaloTrigger/python/L1CaloJets_cff.py new file mode 100644 index 0000000000000..a138f4093e7f4 --- /dev/null +++ b/L1Trigger/L1CaloTrigger/python/L1CaloJets_cff.py @@ -0,0 +1,48 @@ +import FWCore.ParameterSet.Config as cms + +# This cff file is based off of the L1T Phase-2 Menu group using +# the process name "REPR" + + +# -------------------------------------------------------------------------------------------- +# +# ---- Produce the L1EGCrystal clusters using Emulator + +from L1Trigger.L1CaloTrigger.L1EGammaCrystalsEmulatorProducer_cfi import * +L1EGammaClusterEmuProducer.ecalTPEB = cms.InputTag("simEcalEBTriggerPrimitiveDigis","","") + + +# -------------------------------------------------------------------------------------------- +# +# ---- Produce the calibrated tower collection combining Barrel, HGCal, HF + +from L1Trigger.L1CaloTrigger.L1TowerCalibrationProducer_cfi import * +L1TowerCalibrationProducer.L1HgcalTowersInputTag = cms.InputTag("hgcalTowerProducer","HGCalTowerProcessor","") +L1TowerCalibrationProducer.l1CaloTowers = cms.InputTag("L1EGammaClusterEmuProducer","L1CaloTowerCollection","") + + + +# -------------------------------------------------------------------------------------------- +# +# ---- Produce the L1CaloJets + +from L1Trigger.L1CaloTrigger.L1CaloJetProducer_cfi import * +L1CaloJetProducer.l1CaloTowers = cms.InputTag("L1TowerCalibrationProducer","L1CaloTowerCalibratedCollection","") +L1CaloJetProducer.L1CrystalClustersInputTag = cms.InputTag("L1EGammaClusterEmuProducer", "L1EGXtalClusterEmulator","") + + + +# -------------------------------------------------------------------------------------------- +# +# ---- Produce the CaloJet HTT Sums + +from L1Trigger.L1CaloTrigger.L1CaloJetHTTProducer_cfi import * + + + +l1CaloJetsSequence = cms.Sequence( + L1EGammaClusterEmuProducer * + L1TowerCalibrationProducer * + L1CaloJetProducer * + L1CaloJetHTTProducer +) diff --git a/L1Trigger/L1CaloTrigger/test/test_L1CaloJets_cfg.py b/L1Trigger/L1CaloTrigger/test/test_L1CaloJets_cfg.py new file mode 100644 index 0000000000000..498aa095642cc --- /dev/null +++ b/L1Trigger/L1CaloTrigger/test/test_L1CaloJets_cfg.py @@ -0,0 +1,102 @@ +import FWCore.ParameterSet.Config as cms + +from Configuration.StandardSequences.Eras import eras + +process = cms.Process('REPR',eras.Phase2C9) + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') +process.load('SimGeneral.MixingModule.mixNoPU_cfi') +process.load('Configuration.Geometry.GeometryExtended2026D41Reco_cff') +process.load('Configuration.Geometry.GeometryExtended2026D41_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.SimL1Emulator_cff') +process.load('Configuration.StandardSequences.EndOfProcess_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + +process.MessageLogger.categories = cms.untracked.vstring('L1CaloJets', 'FwkReport') +process.MessageLogger.cerr.FwkReport = cms.untracked.PSet( + reportEvery = cms.untracked.int32(1) +) + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(10) ) +#process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) ) + +process.source = cms.Source("PoolSource", + # dasgoclient --query="dataset dataset=/VBFHToTauTau*/Phase2HLTTDR*/FEVT" + #fileNames = cms.untracked.vstring('root://cms-xrd-global.cern.ch//store/mc/PhaseIIMTDTDRAutumn18DR/VBFHToTauTau_M125_14TeV_powheg_pythia8/FEVT/PU200_103X_upgrade2023_realistic_v2-v1/280000/EFC8271A-8026-6A43-AF18-4CB7609B3348.root'), + fileNames = cms.untracked.vstring('root://cms-xrd-global.cern.ch//store/mc/Phase2HLTTDRSummer20ReRECOMiniAOD/VBFHToTauTau_M125_14TeV_powheg_pythia8_correctedGridpack_tuneCP5/FEVT/PU200_111X_mcRun4_realistic_T15_v1-v1/120000/084C8B72-BC64-DE46-801F-D971D5A34F62.root'), + inputCommands = cms.untracked.vstring( + "keep *", + "drop l1tEMTFHitExtras_simEmtfDigis_CSC_HLT", + "drop l1tEMTFHitExtras_simEmtfDigis_RPC_HLT", + "drop l1tEMTFTrackExtras_simEmtfDigis__HLT", + "drop l1tEMTFHit2016Extras_simEmtfDigis_CSC_HLT", + "drop l1tEMTFHit2016Extras_simEmtfDigis_RPC_HLT", + "drop l1tEMTFHit2016s_simEmtfDigis__HLT", + "drop l1tEMTFTrack2016Extras_simEmtfDigis__HLT", + "drop l1tEMTFTrack2016s_simEmtfDigis__HLT", + ) +) + +# All this stuff just runs the various EG algorithms that we are studying + +# ---- Global Tag : +from Configuration.AlCa.GlobalTag import GlobalTag +#process.GlobalTag = GlobalTag(process.GlobalTag, '103X_upgrade2023_realistic_v2', '') +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic', '') + +# Add HCAL Transcoder +process.load('SimCalorimetry.HcalTrigPrimProducers.hcaltpdigi_cff') +process.load('CalibCalorimetry.CaloTPG.CaloTPGTranscoder_cfi') + + +process.L1simulation_step = cms.Path(process.SimL1Emulator) + + + +### Based on: L1Trigger/L1TCommon/test/reprocess_test_10_4_0_mtd5.py +### This code is a portion of what is imported and excludes the 'schedule' portion +### of the two lines below. It makes the test script run! +### from L1Trigger.Configuration.customiseUtils import L1TrackTriggerTracklet +### process = L1TrackTriggerTracklet(process) +process.load('L1Trigger.TrackFindingTracklet.L1HybridEmulationTracks_cff') +process.L1TrackTriggerTracklet_step = cms.Path(process.L1HybridTracksWithAssociators) + + + + +# -------------------------------------------------------------------------------------------- +# +# ---- Load the L1CaloJet sequence designed to accompany process named "REPR" + +process.load('L1Trigger.L1CaloTrigger.L1CaloJets_cff') +process.l1CaloJets = cms.Path(process.l1CaloJetsSequence) + + + +process.Out = cms.OutputModule( "PoolOutputModule", + fileName = cms.untracked.string( "l1caloJetTest.root" ), + fastCloning = cms.untracked.bool( False ), + outputCommands = cms.untracked.vstring( + "keep *_L1EGammaClusterEmuProducer_*_*", + "keep *_L1CaloJetProducer_*_*", + "keep *_L1CaloJetHTTProducer_*_*", + "keep *_TriggerResults_*_*", + "keep *_simHcalTriggerPrimitiveDigis_*_*", + "keep *_EcalEBTrigPrimProducer_*_*", + "keep *_hgcalTowerProducer_*_*" + ) +) + +process.end = cms.EndPath( process.Out ) + + + +#dump_file = open("dump_file.py", "w") +#dump_file.write(process.dumpPython()) + + diff --git a/L1Trigger/Phase2L1ParticleFlow/python/pfClustersFromCombinedCalo_cfi.py b/L1Trigger/Phase2L1ParticleFlow/python/pfClustersFromCombinedCalo_cfi.py index 1e3c3f2bc8612..15cbb2308957f 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/pfClustersFromCombinedCalo_cfi.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/pfClustersFromCombinedCalo_cfi.py @@ -6,7 +6,7 @@ hcalDigis = cms.VInputTag(cms.InputTag('simHcalTriggerPrimitiveDigis')), hcalDigisBarrel = cms.bool(False), hcalDigisHF = cms.bool(True), - phase2barrelCaloTowers = cms.VInputTag(cms.InputTag("L1EGammaClusterEmuProducer",)), + phase2barrelCaloTowers = cms.VInputTag(cms.InputTag("L1EGammaClusterEmuProducer","L1CaloTowerCollection","")), hcalHGCTowers = cms.VInputTag(cms.InputTag("hgcalTowerProducer:HGCalTowerProcessor") ), hcalHGCTowersHadOnly = cms.bool(False), # take also EM part from towers emCorrector = cms.string(""), # no need to correct further