From 071ddb6475b9746736ae6b9a71043901d4fe106e Mon Sep 17 00:00:00 2001 From: Vladimir Date: Thu, 25 Jun 2020 07:30:49 +0200 Subject: [PATCH] Addressed comments of #29904 (b4fdcc3) --- .../interface/ParametricCalibration.h | 28 +- .../L1EGammaCrystalsEmulatorProducer.cc | 598 +++++++++--------- .../plugins/L1EGammaEEProducer.cc | 51 +- .../python/L1TowerCalibrationProducer_cfi.py | 1 - .../python/l1EGammaEEProducer_cfi.py | 2 +- .../L1CaloTrigger/src/L1EGammaEECalibrator.cc | 33 +- .../src/ParametricCalibration.cc | 44 ++ 7 files changed, 376 insertions(+), 381 deletions(-) diff --git a/L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h b/L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h index 9632f5540d972..993c0fedf8230 100644 --- a/L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h +++ b/L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h @@ -2,6 +2,8 @@ #define L1Trigger_L1CaloTrigger_ParametricCalibration_h #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/Utilities/interface/Exception.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" #include #include #include @@ -10,31 +12,9 @@ namespace l1tp2 { class ParametricCalibration { public: ParametricCalibration() {} - ParametricCalibration(const edm::ParameterSet &cpset) { - std::vector etaBins = cpset.getParameter>("etaBins"); - std::vector ptBins = cpset.getParameter>("ptBins"); - std::vector scale = cpset.getParameter>("scale"); - etas.insert(etas.end(), etaBins.begin(), etaBins.end()); - pts.insert(pts.end(), ptBins.begin(), ptBins.end()); - scales.insert(scales.end(), scale.begin(), scale.end()); - if (cpset.existsAs>("ptMin")) { - std::vector ptMin = cpset.getParameter>("ptMin"); - ptMins.insert(ptMins.end(), ptMin.begin(), ptMin.end()); - } else { - float ptMin = cpset.existsAs("ptMin") ? cpset.getParameter("ptMin") : 0; - ptMins = std::vector(etaBins.size(), ptMin); - } - if (cpset.existsAs>("ptMax")) { - std::vector ptMax = cpset.getParameter>("ptMax"); - ptMaxs.insert(ptMaxs.end(), ptMax.begin(), ptMax.end()); - } else { - ptMaxs = std::vector(etaBins.size(), 1e6); - } + ParametricCalibration(const edm::ParameterSet& cpset); + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); - if (pts.size() * etas.size() != scales.size()) - throw cms::Exception("Configuration", - "Bad number of calibration scales, pts.size() * etas.size() != scales.size()"); - } float operator()(const float pt, const float abseta) const { int ptBin = -1; for (unsigned int i = 0, n = pts.size(); i < n; ++i) { diff --git a/L1Trigger/L1CaloTrigger/plugins/L1EGammaCrystalsEmulatorProducer.cc b/L1Trigger/L1CaloTrigger/plugins/L1EGammaCrystalsEmulatorProducer.cc index 88be917079c29..cfc9ed07e285e 100644 --- a/L1Trigger/L1CaloTrigger/plugins/L1EGammaCrystalsEmulatorProducer.cc +++ b/L1Trigger/L1CaloTrigger/plugins/L1EGammaCrystalsEmulatorProducer.cc @@ -54,9 +54,14 @@ Description: Produces crystal clusters using crystal-level information and hardw #include "L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h" #include "L1Trigger/L1TCalorimeter/interface/CaloTools.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" static constexpr bool do_brem = true; +static constexpr int n_eta_bins = 2; +static constexpr int n_borders_phi = 18; +static constexpr int n_borders_eta = 18; +static constexpr int n_clusters_max = 5; static constexpr int n_clusters_link = 3; static constexpr int n_clusters_4link = 4 * 3; static constexpr int n_crystals_towerEta = 5; @@ -70,12 +75,24 @@ static constexpr int n_links_card = 4; static constexpr int n_links_GCTcard = 48; static constexpr float ECAL_eta_range = 1.4841; static constexpr float half_crystal_size = 0.00873; -static constexpr float PiValue = 3.14159; +static constexpr float slideIsoPtThreshold = 80; static constexpr float a0_80 = 0.85, a1_80 = 0.0080, a0 = 0.21; // passes_iso static constexpr float b0 = 0.38, b1 = 1.9, b2 = 0.05; //passes_looseTkiso static constexpr float c0_ss = 0.94, c1_ss = 0.052, c2_ss = 0.044; //passes_ss static constexpr float d0 = 0.96, d1 = 0.0003; //passes_photon static constexpr float e0_looseTkss = 0.944, e1_looseTkss = 0.65, e2_looseTkss = 0.4; //passes_looseTkss +static constexpr float cut_500_MeV = 0.5; + +// absolue IDs range from 0-33 +// 0-16 are iEta -17 to -1 +// 17 to 33 are iEta 1 to 17 +static constexpr int toweriEta_fromAbsoluteID_shift = 16; + +// absolue IDs range from 0-71. +// To align with detector tower IDs (1 - n_towers_Phi) +// shift all indices by 37 and loop over after 72 +static constexpr int toweriPhi_fromAbsoluteID_shift_lowerHalf = 37; +static constexpr int toweriPhi_fromAbsoluteID_shift_upperHalf = 35; float getEta_fromL2LinkCardTowerCrystal(int link, int card, int tower, int crystal) { int etaID = n_crystals_towerEta * (n_towers_cardEta * ((link / n_links_card) % 2) + (tower % n_towers_cardEta)) + @@ -87,8 +104,8 @@ float getEta_fromL2LinkCardTowerCrystal(int link, int card, int tower, int cryst float getPhi_fromL2LinkCardTowerCrystal(int link, int card, int tower, int crystal) { int phiID = n_crystals_towerPhi * ((card * 24) + (n_links_card * (link / 8)) + (tower / n_towers_cardEta)) + crystal / n_crystals_towerPhi; - float size_cell = 2 * PiValue / (n_crystals_towerPhi * n_towers_Phi); - return phiID * size_cell - PiValue + half_crystal_size; + float size_cell = 2 * M_PI / (n_crystals_towerPhi * n_towers_Phi); + return phiID * size_cell - M_PI + half_crystal_size; } int getCrystal_etaID(float eta) { @@ -97,30 +114,6 @@ int getCrystal_etaID(float eta) { return etaID; } -int get_insideEta(float eta, float centereta) { - float size_cell = 2 * ECAL_eta_range / (n_crystals_towerEta * n_towers_Eta); - if (eta - centereta < -1.0 * size_cell / 4) - return 0; - else if (eta - centereta < 0) - return 1; - else if (eta - centereta < 1.0 * size_cell / 4) - return 2; - else - return 3; -} - -int get_insidePhi(float phi, float centerphi) { - float size_cell = 2 * PiValue / (n_crystals_towerPhi * n_towers_Phi); - if (phi - centerphi < -1 * size_cell / 4) - return 0; - else if (phi - centerphi < 0) - return 1; - else if (phi - centerphi < 1 * size_cell / 4) - return 2; - else - return 3; -} - int convert_L2toL1_link(int link) { return link % n_links_card; } int convert_L2toL1_tower(int tower) { return tower; } @@ -128,8 +121,8 @@ int convert_L2toL1_tower(int tower) { return tower; } int convert_L2toL1_card(int card, int link) { return card * n_clusters_4link + link / n_links_card; } int getCrystal_phiID(float phi) { - float size_cell = 2 * PiValue / (n_crystals_towerPhi * n_towers_Phi); - int phiID = int((phi + PiValue) / size_cell); + float size_cell = 2 * M_PI / (n_crystals_towerPhi * n_towers_Phi); + int phiID = int((phi + M_PI) / size_cell); return phiID; } @@ -140,29 +133,23 @@ int getTower_absoluteEtaID(float eta) { } int getTower_absolutePhiID(float phi) { - float size_cell = 2 * PiValue / n_towers_Phi; - int phiID = int((phi + PiValue) / size_cell); + float size_cell = 2 * M_PI / n_towers_Phi; + int phiID = int((phi + M_PI) / size_cell); return phiID; } -// absolue IDs range from 0-33 -// 0-16 are iEta -17 to -1 -// 17 to 33 are iEta 1 to 17 int getToweriEta_fromAbsoluteID(int id) { if (id < n_towers_cardEta) return id - n_towers_cardEta; else - return id - 16; + return id - toweriEta_fromAbsoluteID_shift; } -// absolue IDs range from 0-71. -// To align with detector tower IDs (1 - n_towers_Phi) -// shift all indices by 37 and loop over after 72 int getToweriPhi_fromAbsoluteID(int id) { if (id < n_towers_Phi / 2) - return id + 37; + return id + toweriPhi_fromAbsoluteID_shift_lowerHalf; else - return id - 35; + return id - toweriPhi_fromAbsoluteID_shift_upperHalf; } float getTowerEta_fromAbsoluteID(int id) { @@ -172,8 +159,8 @@ float getTowerEta_fromAbsoluteID(int id) { } float getTowerPhi_fromAbsoluteID(int id) { - float size_cell = 2 * PiValue / n_towers_Phi; - float phi = (id * size_cell) - PiValue + 0.5 * size_cell; + float size_cell = 2 * M_PI / n_towers_Phi; + float phi = (id * size_cell) - M_PI + 0.5 * size_cell; return phi; } @@ -240,7 +227,6 @@ class L1EGCrystalClusterEmulatorProducer : public edm::stream::EDProducer<> { private: void produce(edm::Event&, const edm::EventSetup&) override; - bool passes_he(float pt, float he); bool passes_ss(float pt, float ss); bool passes_photon(float pt, float pss); bool passes_iso(float pt, float iso); @@ -261,60 +247,74 @@ class L1EGCrystalClusterEmulatorProducer : public edm::stream::EDProducer<> { const HcalTopology* hcTopology_; struct mycluster { - float c2x2; - float c2x5; - float c5x5; - int cshowershape; - int cshowershapeloosetk; - float cvalueshowershape; - int cphotonshowershape; - float cpt; // ECAL pt - int cbrem; // if brem corrections were applied - float cWeightedEta; - float cWeightedPhi; - float ciso; // pt of cluster divided by 7x7 ECAL towers - float chovere; // 5x5 HCAL towers divided by the ECAL cluster pt - float craweta; // coordinates between -1.44 and 1.44 - float crawphi; // coordinates between -pi and pi - float chcal; // 5x5 HCAL towers - float ceta; // eta ID in the whole detector (between 0 and 5*34-1) - float cphi; // phi ID in the whole detector (between 0 and 5*72-1) - int ccrystalid; // crystal ID inside tower (between 0 and 24) - int cinsidecrystalid; - int ctowerid; // tower ID inside card (between 0 and 4*n_towers_cardEta-1) + float c2x2_; + float c2x5_; + float c5x5_; + int cshowershape_; + int cshowershapeloosetk_; + float cvalueshowershape_; + int cphotonshowershape_; + float cpt; // ECAL pt + int cbrem_; // if brem corrections were applied + float cWeightedEta_; + float cWeightedPhi_; + float ciso_; // pt of cluster divided by 7x7 ECAL towers + float chovere_; // 5x5 HCAL towers divided by the ECAL cluster pt + float craweta_; // coordinates between -1.44 and 1.44 + float crawphi_; // coordinates between -pi and pi + float chcal_; // 5x5 HCAL towers + float ceta_; // eta ID in the whole detector (between 0 and 5*34-1) + float cphi_; // phi ID in the whole detector (between 0 and 5*72-1) + int ccrystalid_; // crystal ID inside tower (between 0 and 24) + int cinsidecrystalid_; + int ctowerid_; // tower ID inside card (between 0 and 4*n_towers_cardEta-1) }; - bool order_clusters(mycluster c1, mycluster c2) { return c1.cpt < c2.cpt; } - class SimpleCaloHit { - public: - EBDetId id; - HcalDetId id_hcal; - GlobalVector position; // As opposed to GlobalPoint, so we can add them (for weighted average) - float energy = 0.; - bool used = false; - bool stale = false; // Hits become stale once used in clustering algorithm to prevent overlap in clusters - bool isEndcapHit = false; // If using endcap, we won't be using integer crystal indices + private: + float pt_ = 0; + float energy_ = 0.; + bool isEndcapHit_ = false; // If using endcap, we won't be using integer crystal indices + bool stale_ = false; // Hits become stale once used in clustering algorithm to prevent overlap in clusters + bool used_ = false; + GlobalVector position_; // As opposed to GlobalPoint, so we can add them (for weighted average) + HcalDetId id_hcal_; + EBDetId id_; + public: // tool functions - inline float pt() const { return (position.mag2() > 0) ? energy * sin(position.theta()) : 0.; }; - inline float deta(SimpleCaloHit& other) const { return position.eta() - other.position.eta(); }; + inline void setPt() { pt_ = (position_.mag2() > 0) ? energy_ * sin(position_.theta()) : 0; }; + inline void setEnergy(float et) { energy_ = et / sin(position_.theta()); }; + inline void setIsEndcapHit(bool isEC) { isEndcapHit_ = isEC; }; + inline void setUsed(bool isUsed) { used_ = isUsed; }; + inline void setPosition(const GlobalVector& pos) { position_ = pos; }; + inline void setIdHcal(const HcalDetId& idhcal) { id_hcal_ = idhcal; }; + inline void setId(const EBDetId& id) { id_ = id; }; + + inline float pt() const { return pt_; }; + inline float energy() const { return energy_; }; + inline bool isEndcapHit() const { return isEndcapHit_; }; + inline bool used() const { return used_; }; + inline const GlobalVector& position() const { return position_; }; + inline const EBDetId& id() const { return id_; }; + + inline float deta(SimpleCaloHit& other) const { return position_.eta() - other.position().eta(); }; int dieta(SimpleCaloHit& other) const { - if (isEndcapHit || other.isEndcapHit) + if (isEndcapHit_ || other.isEndcapHit()) return 9999; // We shouldn't compare integer indices in endcap, the map is not linear - if (id.ieta() * other.id.ieta() > 0) - return id.ieta() - other.id.ieta(); - return id.ieta() - other.id.ieta() - 1; + if (id_.ieta() * other.id().ieta() > 0) + return id_.ieta() - other.id().ieta(); + return id_.ieta() - other.id().ieta() - 1; }; inline float dphi(SimpleCaloHit& other) const { - return reco::deltaPhi(static_cast(position.phi()), static_cast(other.position.phi())); + return reco::deltaPhi(static_cast(position_.phi()), static_cast(other.position().phi())); }; int diphi(SimpleCaloHit& other) const { - if (isEndcapHit || other.isEndcapHit) + if (isEndcapHit_ || other.isEndcapHit()) return 9999; // We shouldn't compare integer indices in endcap, the map is not linear // Logic from EBDetId::distancePhi() without the abs() - int PI = 180; - int result = id.iphi() - other.id.iphi(); + static constexpr int PI = 180; + int result = id().iphi() - other.id().iphi(); while (result > PI) result -= 2 * PI; while (result <= -PI) @@ -324,13 +324,11 @@ class L1EGCrystalClusterEmulatorProducer : public edm::stream::EDProducer<> { float distanceTo(SimpleCaloHit& other) const { // Treat position as a point, measure 3D distance // This is used for endcap hits, where we don't have a rectangular mapping - return (position - other.position).mag(); + return (position() - other.position()).mag(); }; bool operator==(SimpleCaloHit& other) const { - if (id == other.id && position == other.position && energy == other.energy && isEndcapHit == other.isEndcapHit) - return true; - - return false; + return (id_ == other.id() && position() == other.position() && energy_ == other.energy() && + isEndcapHit_ == other.isEndcapHit()); }; }; }; @@ -369,20 +367,21 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: iEvent.getByToken(ecalTPEBToken_, pcalohits); std::vector ecalhits; - for (auto& hit : *pcalohits.product()) { + for (const auto& hit : *pcalohits.product()) { if (hit.encodedEt() > 0) // hit.encodedEt() returns an int corresponding to 2x the crystal Et { - float et = hit.encodedEt() / - 8.; // Et is 10 bit, by keeping the ADC saturation Et at 120 GeV it means that you have to divide by 8 - if (et < 0.5) + // Et is 10 bit, by keeping the ADC saturation Et at 120 GeV it means that you have to divide by 8 + float et = hit.encodedEt() / 8.; + if (et < cut_500_MeV) continue; // keep the 500 MeV ET Cut auto cell = ebGeometry->getGeometry(hit.id()); SimpleCaloHit ehit; - ehit.id = hit.id(); - ehit.position = GlobalVector(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z()); - ehit.energy = et / sin(ehit.position.theta()); + ehit.setId(hit.id()); + ehit.setPosition(GlobalVector(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z())); + ehit.setEnergy(et); + ehit.setPt(); ecalhits.push_back(ehit); } } @@ -391,24 +390,27 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: std::vector hcalhits; edm::Handle > hbhecoll; iEvent.getByToken(hcalTPToken_, hbhecoll); - for (auto& hit : *hbhecoll.product()) { - //if ( hit.SOI_compressedEt() == 0 ) continue; + for (const auto& hit : *hbhecoll.product()) { float et = decoder_->hcaletValue(hit.id(), hit.t0()); if (et <= 0) continue; if (!(hcTopology_->validHT(hit.id()))) { - std::cout << " -- Hcal hit DetID not present in HCAL Geom: " << hit.id() << std::endl; + LogError("L1EGCrystalClusterEmulatorProducer") + << " -- Hcal hit DetID not present in HCAL Geom: " << hit.id() << std::endl; + throw cms::Exception("L1EGCrystalClusterEmulatorProducer"); continue; } std::vector hcId = theTrigTowerGeometry.detIds(hit.id()); if (hcId.empty()) { - std::cout << "Cannot find any HCalDetId corresponding to " << hit.id() << std::endl; + LogError("L1EGCrystalClusterEmulatorProducer") + << "Cannot find any HCalDetId corresponding to " << hit.id() << std::endl; + throw cms::Exception("L1EGCrystalClusterEmulatorProducer"); continue; } if (hcId[0].subdetId() > 1) continue; GlobalVector hcal_tp_position = GlobalVector(0., 0., 0.); - for (auto& hcId_i : hcId) { + for (const auto& hcId_i : hcId) { if (hcId_i.subdetId() > 1) continue; auto cell = hbGeometry->getGeometry(hcId_i); @@ -419,11 +421,11 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: break; } SimpleCaloHit hhit; - hhit.id = hit.id(); - hhit.id_hcal = hit.id(); - hhit.position = hcal_tp_position; - //float et = hit.SOI_compressedEt() / 2.; - hhit.energy = et / sin(hhit.position.theta()); + hhit.setId(hit.id()); + hhit.setIdHcal(hit.id()); + hhit.setPosition(hcal_tp_position); + hhit.setEnergy(et); + hhit.setPt(); hcalhits.push_back(hhit); } @@ -432,28 +434,20 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: //******************************************************************* // Definition of L1 outputs - float ECAL_tower_L1Card[n_links_card][n_towers_cardEta] - [n_towers_halfPhi]; // 36 L1 cards send each 4 links with 17 towers - float HCAL_tower_L1Card[n_links_card][n_towers_cardEta] - [n_towers_halfPhi]; // 36 L1 cards send each 4 links with 17 towers - int iEta_tower_L1Card[n_links_card][n_towers_cardEta] - [n_towers_halfPhi]; // 36 L1 cards send each 4 links with 17 towers - int iPhi_tower_L1Card[n_links_card][n_towers_cardEta] - [n_towers_halfPhi]; // 36 L1 cards send each 4 links with 17 towers - float energy_cluster_L1Card[n_links_card][n_clusters_link] - [n_towers_halfPhi]; // 36 L1 cards send each 4 links with 3 clusters - int brem_cluster_L1Card[n_links_card][n_clusters_link] - [n_towers_halfPhi]; // 36 L1 cards send each 4 links with 3 clusters - int towerID_cluster_L1Card[n_links_card][n_clusters_link] - [n_towers_halfPhi]; // 36 L1 cards send each 4 links with 3 clusters - int crystalID_cluster_L1Card[n_links_card][n_clusters_link] - [n_towers_halfPhi]; // 36 L1 cards send each 4 links with 3 clusters - int showerShape_cluster_L1Card[n_links_card][n_clusters_link] - [n_towers_halfPhi]; // 36 L1 cards send each 4 links with 3 clusters - int showerShapeLooseTk_cluster_L1Card[n_links_card][n_clusters_link] - [n_towers_halfPhi]; // 36 L1 cards send each 4 links with 3 clusters - int photonShowerShape_cluster_L1Card[n_links_card][n_clusters_link] - [n_towers_halfPhi]; // 36 L1 cards send each 4 links with 3 clusters + // 36 L1 cards send each 4 links with 17 towers + float ECAL_tower_L1Card[n_links_card][n_towers_cardEta][n_towers_halfPhi]; + float HCAL_tower_L1Card[n_links_card][n_towers_cardEta][n_towers_halfPhi]; + int iEta_tower_L1Card[n_links_card][n_towers_cardEta][n_towers_halfPhi]; + int iPhi_tower_L1Card[n_links_card][n_towers_cardEta][n_towers_halfPhi]; + // 36 L1 cards send each 4 links with 3 clusters + float energy_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi]; + // 36 L1 cards send each 4 links with 3 clusters + int brem_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi]; + int towerID_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi]; + int crystalID_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi]; + int showerShape_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi]; + int showerShapeLooseTk_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi]; + int photonShowerShape_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi]; for (int ii = 0; ii < n_links_card; ++ii) { for (int jj = 0; jj < n_towers_cardEta; ++jj) { @@ -476,43 +470,46 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: } } - vector - cluster_list[n_towers_halfPhi]; // There is one list of clusters per card. We take the 12 highest pt per card - vector - cluster_list_merged[n_towers_halfPhi]; // After merging the clusters in different regions of a single L1 card + // There is one list of clusters per card. We take the 12 highest pt per card + vector cluster_list[n_towers_halfPhi]; + // After merging the clusters in different regions of a single L1 card + vector cluster_list_merged[n_towers_halfPhi]; for (int cc = 0; cc < n_towers_halfPhi; ++cc) { // Loop over 36 L1 cards - - for (int nregion = 0; nregion < 6; ++nregion) { // Loop over 3x4 etaxphi regions to search for max 5 clusters + // Loop over 3x4 etaxphi regions to search for max 5 clusters + for (int nregion = 0; nregion <= n_clusters_max; ++nregion) { int nclusters = 0; bool build_cluster = true; - while (nclusters < 5 && build_cluster) { // Continue until 5 clusters have been built or there is no cluster left + // Continue until 5 clusters have been built or there is no cluster left + while (nclusters < n_clusters_max && build_cluster) { build_cluster = false; SimpleCaloHit centerhit; for (const auto& hit : ecalhits) { - if (getCrystal_phiID(hit.position.phi()) <= getPhiMax_card(cc) && - getCrystal_phiID(hit.position.phi()) >= getPhiMin_card(cc) && - getCrystal_etaID(hit.position.eta()) <= getEtaMax_card(cc) && - getCrystal_etaID(hit.position.eta()) >= getEtaMin_card(cc)) { // Check that the hit is in the good card - if (getCrystal_etaID(hit.position.eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) && - getCrystal_etaID(hit.position.eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion && - !hit.used && hit.pt() >= 1.0 && hit.pt() > centerhit.pt()) // 3 towers x 5 crystals - { // Highest hit in good region with pt>1 and not used in any other cluster - centerhit = hit; - build_cluster = true; - } + if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) && + getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) && + getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) && + getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && + // Check that the hit is in the good card + getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) && + getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion && + !hit.used() && hit.pt() >= 1.0 && hit.pt() > centerhit.pt()) // 3 towers x 5 crystals + { + // Highest hit in good region with pt>1 and not used in any other cluster + centerhit = hit; + build_cluster = true; } } if (build_cluster) nclusters++; - if (build_cluster && nclusters > 0 && nclusters < 6) { // Use only the 5 most energetic clusters + // Use only the 5 most energetic clusters + if (build_cluster && nclusters > 0 && nclusters <= n_clusters_max) { mycluster mc1; mc1.cpt = 0.0; - mc1.cWeightedEta = 0.0; - mc1.cWeightedPhi = 0.0; + mc1.cWeightedEta_ = 0.0; + mc1.cWeightedPhi_ = 0.0; float leftlobe = 0; float rightlobe = 0; float e5x5 = 0; @@ -530,12 +527,12 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: float e2x2_4 = 0; float n2x2_4 = 0; for (auto& hit : ecalhits) { - if (getCrystal_phiID(hit.position.phi()) <= getPhiMax_card(cc) && - getCrystal_phiID(hit.position.phi()) >= getPhiMin_card(cc) && - getCrystal_etaID(hit.position.eta()) <= getEtaMax_card(cc) && - getCrystal_etaID(hit.position.eta()) >= getEtaMin_card(cc) && hit.pt() > 0 && - getCrystal_etaID(hit.position.eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) && - getCrystal_etaID(hit.position.eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion) { + if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) && + getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) && + getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) && + getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && hit.pt() > 0 && + getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) && + getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion) { if (abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) > 2 && hit.diphi(centerhit) <= 7) { rightlobe += hit.pt(); } @@ -543,95 +540,91 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: leftlobe += hit.pt(); } if (abs(hit.dieta(centerhit)) <= 2 && abs(hit.diphi(centerhit)) <= 2) { - e5x5 += hit.energy; + e5x5 += hit.energy(); n5x5++; } if ((hit.dieta(centerhit) == 1 or hit.dieta(centerhit) == 0) && (hit.diphi(centerhit) == 1 or hit.diphi(centerhit) == 0)) { - e2x2_1 += hit.energy; + e2x2_1 += hit.energy(); n2x2_1++; } if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) && (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == 1)) { - e2x2_2 += hit.energy; + e2x2_2 += hit.energy(); n2x2_2++; } if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == 1) && (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == -1)) { - e2x2_3 += hit.energy; + e2x2_3 += hit.energy(); n2x2_3++; } if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) && (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == -1)) { - e2x2_4 += hit.energy; + e2x2_4 += hit.energy(); n2x2_4++; } if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == 1) && abs(hit.diphi(centerhit)) <= 2) { - e2x5_1 += hit.energy; + e2x5_1 += hit.energy(); n2x5_1++; } if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) && abs(hit.diphi(centerhit)) <= 2) { - e2x5_2 += hit.energy; + e2x5_2 += hit.energy(); n2x5_2++; } } - if (getCrystal_phiID(hit.position.phi()) <= getPhiMax_card(cc) && - getCrystal_phiID(hit.position.phi()) >= getPhiMin_card(cc) && - getCrystal_etaID(hit.position.eta()) <= getEtaMax_card(cc) && - getCrystal_etaID(hit.position.eta()) >= getEtaMin_card(cc) && !hit.used && hit.pt() > 0 && + if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) && + getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) && + getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) && + getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && !hit.used() && hit.pt() > 0 && abs(hit.dieta(centerhit)) <= 1 && abs(hit.diphi(centerhit)) <= 2 && - getCrystal_etaID(hit.position.eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) && - getCrystal_etaID(hit.position.eta()) >= - getEtaMin_card(cc) + - n_crystals_3towers * - nregion) { // clusters 3x5 in etaxphi using only the hits in the corresponding card and in the corresponding 3x4 region - hit.used = true; + getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) && + getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion) { + // clusters 3x5 in etaxphi using only the hits in the corresponding card and in the corresponding 3x4 region + hit.setUsed(true); mc1.cpt += hit.pt(); - mc1.cWeightedEta += float(hit.pt()) * float(hit.position.eta()); - mc1.cWeightedPhi = mc1.cWeightedPhi + (float(hit.pt()) * float(hit.position.phi())); + mc1.cWeightedEta_ += float(hit.pt()) * float(hit.position().eta()); + mc1.cWeightedPhi_ = mc1.cWeightedPhi_ + (float(hit.pt()) * float(hit.position().phi())); } } if (do_brem && (rightlobe > 0.10 * mc1.cpt or leftlobe > 0.10 * mc1.cpt)) { for (auto& hit : ecalhits) { - if (getCrystal_phiID(hit.position.phi()) <= getPhiMax_card(cc) && - getCrystal_phiID(hit.position.phi()) >= getPhiMin_card(cc) && - getCrystal_etaID(hit.position.eta()) <= getEtaMax_card(cc) && - getCrystal_etaID(hit.position.eta()) >= getEtaMin_card(cc) && hit.pt() > 0 && - getCrystal_etaID(hit.position.eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) && - getCrystal_etaID(hit.position.eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion && - !hit.used) { + if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) && + getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) && + getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) && + getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && hit.pt() > 0 && + getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) && + getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion && + !hit.used()) { if (rightlobe > 0.10 * mc1.cpt && (leftlobe < 0.10 * mc1.cpt or rightlobe > leftlobe) && abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) > 2 && hit.diphi(centerhit) <= 7) { mc1.cpt += hit.pt(); - hit.used = true; - mc1.cbrem = 1; + hit.setUsed(true); + mc1.cbrem_ = 1; } if (leftlobe > 0.10 * mc1.cpt && (rightlobe < 0.10 * mc1.cpt or leftlobe >= rightlobe) && abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) < -2 && hit.diphi(centerhit) >= -7) { mc1.cpt += hit.pt(); - hit.used = true; - mc1.cbrem = 1; + hit.setUsed(true); + mc1.cbrem_ = 1; } } } } - mc1.c5x5 = e5x5; - mc1.c2x5 = e2x5_1; - if (e2x5_2 > mc1.c2x5) - mc1.c2x5 = e2x5_2; - mc1.c2x2 = e2x2_1; - if (e2x2_2 > mc1.c2x2) - mc1.c2x2 = e2x2_2; - if (e2x2_3 > mc1.c2x2) - mc1.c2x2 = e2x2_3; - if (e2x2_4 > mc1.c2x2) - mc1.c2x2 = e2x2_4; - mc1.cWeightedEta = mc1.cWeightedEta / mc1.cpt; - mc1.cWeightedPhi = mc1.cWeightedPhi / mc1.cpt; - mc1.ceta = getCrystal_etaID(centerhit.position.eta()); - mc1.cphi = getCrystal_phiID(centerhit.position.phi()); - mc1.crawphi = centerhit.position.phi(); - mc1.craweta = centerhit.position.eta(); + mc1.c5x5_ = e5x5; + mc1.c2x5_ = max(e2x5_1, e2x5_2); + mc1.c2x2_ = e2x2_1; + if (e2x2_2 > mc1.c2x2_) + mc1.c2x2_ = e2x2_2; + if (e2x2_3 > mc1.c2x2_) + mc1.c2x2_ = e2x2_3; + if (e2x2_4 > mc1.c2x2_) + mc1.c2x2_ = e2x2_4; + mc1.cWeightedEta_ = mc1.cWeightedEta_ / mc1.cpt; + mc1.cWeightedPhi_ = mc1.cWeightedPhi_ / mc1.cpt; + mc1.ceta_ = getCrystal_etaID(centerhit.position().eta()); + mc1.cphi_ = getCrystal_phiID(centerhit.position().phi()); + mc1.crawphi_ = centerhit.position().phi(); + mc1.craweta_ = centerhit.position().eta(); cluster_list[cc].push_back(mc1); } // End if 5 clusters per region } // End while to find the 5 clusters @@ -641,24 +634,24 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: // Merge clusters from different regions for (unsigned int jj = 0; jj < unsigned(cluster_list[cc].size()); ++jj) { for (unsigned int kk = jj + 1; kk < unsigned(cluster_list[cc].size()); ++kk) { - if (std::abs(cluster_list[cc][jj].ceta - cluster_list[cc][kk].ceta) < 2 && - std::abs(cluster_list[cc][jj].cphi - cluster_list[cc][kk].cphi) < 2) { //Diagonale + exact neighbors + if (std::abs(cluster_list[cc][jj].ceta_ - cluster_list[cc][kk].ceta_) < 2 && + std::abs(cluster_list[cc][jj].cphi_ - cluster_list[cc][kk].cphi_) < 2) { //Diagonale + exact neighbors if (cluster_list[cc][kk].cpt > cluster_list[cc][jj].cpt) { cluster_list[cc][kk].cpt += cluster_list[cc][jj].cpt; - cluster_list[cc][kk].c5x5 += cluster_list[cc][jj].c5x5; - cluster_list[cc][kk].c2x5 += cluster_list[cc][jj].c2x5; + cluster_list[cc][kk].c5x5_ += cluster_list[cc][jj].c5x5_; + cluster_list[cc][kk].c2x5_ += cluster_list[cc][jj].c2x5_; cluster_list[cc][jj].cpt = 0; - cluster_list[cc][jj].c5x5 = 0; - cluster_list[cc][jj].c2x5 = 0; - cluster_list[cc][jj].c2x2 = 0; + cluster_list[cc][jj].c5x5_ = 0; + cluster_list[cc][jj].c2x5_ = 0; + cluster_list[cc][jj].c2x2_ = 0; } else { cluster_list[cc][jj].cpt += cluster_list[cc][kk].cpt; - cluster_list[cc][jj].c5x5 += cluster_list[cc][kk].c5x5; - cluster_list[cc][jj].c2x5 += cluster_list[cc][kk].c2x5; + cluster_list[cc][jj].c5x5_ += cluster_list[cc][kk].c5x5_; + cluster_list[cc][jj].c2x5_ += cluster_list[cc][kk].c2x5_; cluster_list[cc][kk].cpt = 0; - cluster_list[cc][kk].c2x2 = 0; - cluster_list[cc][kk].c2x5 = 0; - cluster_list[cc][kk].c5x5 = 0; + cluster_list[cc][kk].c2x2_ = 0; + cluster_list[cc][kk].c2x5_ = 0; + cluster_list[cc][kk].c5x5_ = 0; } } } @@ -666,7 +659,7 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: cluster_list[cc][jj].cpt = cluster_list[cc][jj].cpt * calib_(cluster_list[cc][jj].cpt, - std::abs(cluster_list[cc][jj].craweta)); //Mark's calibration as a function of eta and pt + std::abs(cluster_list[cc][jj].craweta_)); //Mark's calibration as a function of eta and pt cluster_list_merged[cc].push_back(cluster_list[cc][jj]); } } @@ -677,23 +670,23 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: // Fill cluster information in the arrays. We keep max 12 clusters (distributed between 4 links) for (unsigned int jj = 0; jj < unsigned(cluster_list_merged[cc].size()) && jj < n_clusters_4link; ++jj) { crystalID_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = - getCrystalIDInTower(cluster_list_merged[cc][jj].ceta, cluster_list_merged[cc][jj].cphi); + getCrystalIDInTower(cluster_list_merged[cc][jj].ceta_, cluster_list_merged[cc][jj].cphi_); towerID_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = - getTowerID(cluster_list_merged[cc][jj].ceta, cluster_list_merged[cc][jj].cphi); + getTowerID(cluster_list_merged[cc][jj].ceta_, cluster_list_merged[cc][jj].cphi_); energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = cluster_list_merged[cc][jj].cpt; - brem_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = cluster_list_merged[cc][jj].cbrem; + brem_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = cluster_list_merged[cc][jj].cbrem_; if (passes_ss(cluster_list_merged[cc][jj].cpt, - cluster_list_merged[cc][jj].c2x5 / cluster_list_merged[cc][jj].c5x5)) + cluster_list_merged[cc][jj].c2x5_ / cluster_list_merged[cc][jj].c5x5_)) showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1; else showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0; if (passes_looseTkss(cluster_list_merged[cc][jj].cpt, - cluster_list_merged[cc][jj].c2x5 / cluster_list_merged[cc][jj].c5x5)) + cluster_list_merged[cc][jj].c2x5_ / cluster_list_merged[cc][jj].c5x5_)) showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1; else showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0; if (passes_photon(cluster_list_merged[cc][jj].cpt, - cluster_list_merged[cc][jj].c2x2 / cluster_list_merged[cc][jj].c2x5)) + cluster_list_merged[cc][jj].c2x2_ / cluster_list_merged[cc][jj].c2x5_)) photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1; else photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0; @@ -701,19 +694,19 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: // Loop over calo ecal hits to get the ECAL towers. Take only hits that have not been used to make clusters for (const auto& hit : ecalhits) { - if (getCrystal_phiID(hit.position.phi()) <= getPhiMax_card(cc) && - getCrystal_phiID(hit.position.phi()) >= getPhiMin_card(cc) && - getCrystal_etaID(hit.position.eta()) <= getEtaMax_card(cc) && - getCrystal_etaID(hit.position.eta()) >= getEtaMin_card(cc) && - !hit.used) { // Take all the hits inside the card that have not been used yet + if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) && + getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) && + getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) && + getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && + !hit.used()) { // Take all the hits inside the card that have not been used yet for (int jj = 0; jj < n_links_card; ++jj) { // loop over 4 links per card - if ((getCrystal_phiID(hit.position.phi()) / n_crystals_towerPhi) % 4 == jj) { // Go to ID tower modulo 4 + if ((getCrystal_phiID(hit.position().phi()) / n_crystals_towerPhi) % 4 == jj) { // Go to ID tower modulo 4 for (int ii = 0; ii < n_towers_cardEta; ++ii) { //Apply Mark's calibration at the same time (row of the lowest pT, as a function of eta) - if ((getCrystal_etaID(hit.position.eta()) / n_crystals_towerEta) % n_towers_cardEta == ii) { - ECAL_tower_L1Card[jj][ii][cc] += hit.pt() * calib_(0, std::abs(hit.position.eta())); - iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(hit.position.eta()); //hit.id.ieta(); - iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(hit.position.phi()); //hit.id.iphi(); + if ((getCrystal_etaID(hit.position().eta()) / n_crystals_towerEta) % n_towers_cardEta == ii) { + ECAL_tower_L1Card[jj][ii][cc] += hit.pt() * calib_(0, std::abs(hit.position().eta())); + iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(hit.position().eta()); //hit.id().ieta(); + iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(hit.position().phi()); //hit.id().iphi(); } } // end of loop over eta towers } @@ -722,7 +715,7 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: static constexpr float tower_width = 0.0873; for (int jj = 0; jj < n_links_card; ++jj) { for (int ii = 0; ii < n_towers_cardEta; ++ii) { - float phi = getPhiMin_card(cc) * tower_width / n_crystals_towerPhi - PiValue + (jj + 0.5) * tower_width; + float phi = getPhiMin_card(cc) * tower_width / n_crystals_towerPhi - M_PI + (jj + 0.5) * tower_width; float eta = getEtaMin_card(cc) * tower_width / n_crystals_towerEta - n_towers_cardEta * tower_width + (ii + 0.5) * tower_width; iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(eta); @@ -735,17 +728,17 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: // Loop over hcal hits to get the HCAL towers. for (const auto& hit : hcalhits) { - if (getCrystal_phiID(hit.position.phi()) <= getPhiMax_card(cc) && - getCrystal_phiID(hit.position.phi()) >= getPhiMin_card(cc) && - getCrystal_etaID(hit.position.eta()) <= getEtaMax_card(cc) && - getCrystal_etaID(hit.position.eta()) >= getEtaMin_card(cc) && hit.pt() > 0) { + if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) && + getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) && + getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) && + getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && hit.pt() > 0) { for (int jj = 0; jj < n_links_card; ++jj) { - if ((getCrystal_phiID(hit.position.phi()) / n_crystals_towerPhi) % n_links_card == jj) { + if ((getCrystal_phiID(hit.position().phi()) / n_crystals_towerPhi) % n_links_card == jj) { for (int ii = 0; ii < n_towers_cardEta; ++ii) { - if ((getCrystal_etaID(hit.position.eta()) / n_crystals_towerEta) % n_towers_cardEta == ii) { + if ((getCrystal_etaID(hit.position().eta()) / n_crystals_towerEta) % n_towers_cardEta == ii) { HCAL_tower_L1Card[jj][ii][cc] += hit.pt(); - iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(hit.position.eta()); //hit.id.ieta(); - iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(hit.position.phi()); //hit.id.iphi(); + iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(hit.position().eta()); //hit.id().ieta(); + iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(hit.position().phi()); //hit.id().iphi(); } } // end of loop over eta towers } @@ -756,8 +749,8 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: // Give back energy of not used clusters to the towers (if there are more than 12 clusters) for (unsigned int kk = n_clusters_4link; kk < cluster_list_merged[cc].size(); ++kk) { if (cluster_list_merged[cc][kk].cpt > 0) { - ECAL_tower_L1Card[getTower_phiID(cluster_list_merged[cc][kk].cphi)] - [getTower_etaID(cluster_list_merged[cc][kk].ceta)][cc] += cluster_list_merged[cc][kk].cpt; + ECAL_tower_L1Card[getTower_phiID(cluster_list_merged[cc][kk].cphi_)] + [getTower_etaID(cluster_list_merged[cc][kk].ceta_)][cc] += cluster_list_merged[cc][kk].cpt; } } } //End of loop over cards @@ -807,12 +800,12 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: } } - vector cluster_list_L2 - [n_towers_halfPhi]; // There is one list of clusters per equivalent of L1 card. We take the 8 highest pt. + // There is one list of clusters per equivalent of L1 card. We take the 8 highest pt. + vector cluster_list_L2[n_towers_halfPhi]; // Merge clusters on the phi edges - for (int ii = 0; ii < 18; ++ii) { // 18 borders in phi - for (int jj = 0; jj < 2; ++jj) { // 2 eta bins + for (int ii = 0; ii < n_borders_phi; ++ii) { // 18 borders in phi + for (int jj = 0; jj < n_eta_bins; ++jj) { // 2 eta bins int card_left = 2 * ii + jj; int card_right = 2 * ii + jj + 2; if (card_right > 35) @@ -849,30 +842,30 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: // Bremsstrahlung corrections. Merge clusters on the phi edges depending on pT (pt more than 10 percent, dphi leq 5, deta leq 1) if (do_brem) { - for (int ii = 0; ii < 18; ++ii) { // 18 borders in phi - for (int jj = 0; jj < 2; ++jj) { // 2 eta bins + for (int ii = 0; ii < n_borders_phi; ++ii) { // 18 borders in phi + for (int jj = 0; jj < n_eta_bins; ++jj) { // 2 eta bins int card_left = 2 * ii + jj; int card_right = 2 * ii + jj + 2; if (card_right > 35) card_right = card_right - n_towers_halfPhi; - for (int kk = 0; kk < n_clusters_4link; - ++kk) { // 12 clusters in the first card. We check the right side crystalID_cluster_L1Card[kk%4][kk/4][card_left] + // 12 clusters in the first card. We check the right side crystalID_cluster_L1Card[kk%4][kk/4][card_left] + for (int kk = 0; kk < n_clusters_4link; ++kk) { + // if the tower is at the edge there might be brem corrections, whatever the crystal ID if (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 50 && - energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > - 0) { // if the tower is at the edge there might be brem corrections, whatever the crystal ID + energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 0) { for (int ll = 0; ll < n_clusters_4link; ++ll) { // We check the 12 clusters in the card on the right + //Distance of 1 max in eta if (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < n_towers_cardEta && std::abs(5 * (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right]) % n_towers_cardEta + crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] % 5 - 5 * (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left]) % n_towers_cardEta - - crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] % 5) <= - 1) { //Distance of 1 max in eta + crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] % 5) <= 1) { + //Distance of 5 max in phi if (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < n_towers_cardEta && std::abs(5 + crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] / 5 - - (crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] / 5)) <= - 5) { //Distance of 5 max in phi + (crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] / 5)) <= 5) { if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] && energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] > @@ -881,7 +874,8 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right]; energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] = 0; brem_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] = 1; - } // The least energetic cluster is merged to the most energetic one, if its pt is at least ten percent + } + // The least energetic cluster is merged to the most energetic one, if its pt is at least ten percent else if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_right] >= energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_left] && energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_left] > @@ -901,7 +895,7 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: } // Merge clusters on the eta edges - for (int ii = 0; ii < 18; ++ii) { // 18 borders in eta + for (int ii = 0; ii < n_borders_eta; ++ii) { // 18 borders in eta int card_bottom = 2 * ii; int card_top = 2 * ii + 1; for (int kk = 0; kk < n_clusters_4link; ++kk) { // 12 clusters in the first card. We check the top side @@ -937,12 +931,12 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: if (energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii] > 0) { mycluster mc1; mc1.cpt = energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii]; - mc1.cbrem = brem_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii]; - mc1.ctowerid = towerID_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii]; - mc1.ccrystalid = crystalID_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii]; - mc1.cshowershape = showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii]; - mc1.cshowershapeloosetk = showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii]; - mc1.cphotonshowershape = photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii]; + mc1.cbrem_ = brem_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii]; + mc1.ctowerid_ = towerID_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii]; + mc1.ccrystalid_ = crystalID_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii]; + mc1.cshowershape_ = showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii]; + mc1.cshowershapeloosetk_ = showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii]; + mc1.cphotonshowershape_ = photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii]; cluster_list_L2[ii].push_back(mc1); } } @@ -954,11 +948,11 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: for (int ii = 0; ii < n_towers_halfPhi; ++ii) { for (unsigned int jj = 8; jj < n_clusters_4link && jj < cluster_list_L2[ii].size(); ++jj) { if (cluster_list_L2[ii][jj].cpt > 0) { - ECAL_tower_L1Card[cluster_list_L2[ii][jj].ctowerid / n_towers_cardEta] - [cluster_list_L2[ii][jj].ctowerid % n_towers_cardEta][ii] += cluster_list_L2[ii][jj].cpt; + ECAL_tower_L1Card[cluster_list_L2[ii][jj].ctowerid_ / n_towers_cardEta] + [cluster_list_L2[ii][jj].ctowerid_ % n_towers_cardEta][ii] += cluster_list_L2[ii][jj].cpt; cluster_list_L2[ii][jj].cpt = 0; - cluster_list_L2[ii][jj].ctowerid = 0; - cluster_list_L2[ii][jj].ccrystalid = 0; + cluster_list_L2[ii][jj].ctowerid_ = 0; + cluster_list_L2[ii][jj].ccrystalid_ = 0; } } } @@ -966,16 +960,16 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: // Compute isolation (7*7 ECAL towers) and HCAL energy (5x5 HCAL towers) for (int ii = 0; ii < n_towers_halfPhi; ++ii) { // Loop over the new cluster list (stored in 36x8 format) for (unsigned int jj = 0; jj < 8 && jj < cluster_list_L2[ii].size(); ++jj) { - int cluster_etaOfTower_fullDetector = get_towerEta_fromCardTowerInCard(ii, cluster_list_L2[ii][jj].ctowerid); - int cluster_phiOfTower_fullDetector = get_towerPhi_fromCardTowerInCard(ii, cluster_list_L2[ii][jj].ctowerid); + int cluster_etaOfTower_fullDetector = get_towerEta_fromCardTowerInCard(ii, cluster_list_L2[ii][jj].ctowerid_); + int cluster_phiOfTower_fullDetector = get_towerPhi_fromCardTowerInCard(ii, cluster_list_L2[ii][jj].ctowerid_); float hcal_nrj = 0.0; float isolation = 0.0; int ntowers = 0; for (int iii = 0; iii < n_towers_halfPhi; ++iii) { // The clusters have to be added to the isolation for (unsigned int jjj = 0; jjj < 8 && jjj < cluster_list_L2[iii].size(); ++jjj) { if (!(iii == ii && jjj == jj)) { - int cluster2_eta = get_towerEta_fromCardTowerInCard(iii, cluster_list_L2[iii][jjj].ctowerid); - int cluster2_phi = get_towerPhi_fromCardTowerInCard(iii, cluster_list_L2[iii][jjj].ctowerid); + int cluster2_eta = get_towerEta_fromCardTowerInCard(iii, cluster_list_L2[iii][jjj].ctowerid_); + int cluster2_phi = get_towerPhi_fromCardTowerInCard(iii, cluster_list_L2[iii][jjj].ctowerid_); if (abs(cluster2_eta - cluster_etaOfTower_fullDetector) <= 2 && (abs(cluster2_phi - cluster_phiOfTower_fullDetector) <= 2 or abs(cluster2_phi - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) { @@ -985,38 +979,38 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: } } for (int kk = 0; kk < n_towers_halfPhi; ++kk) { // 36 cards - for (int ll = 0; ll < 4; ++ll) { // 4 links per card + for (int ll = 0; ll < n_links_card; ++ll) { // 4 links per card for (int mm = 0; mm < n_towers_cardEta; ++mm) { // 17 towers per link int etaOftower_fullDetector = get_towerEta_fromCardLinkTower(kk, ll, mm); int phiOftower_fullDetector = get_towerPhi_fromCardLinkTower(kk, ll, mm); // First do ECAL + // The towers are within 3. Needs to stitch the two phi sides together if (abs(etaOftower_fullDetector - cluster_etaOfTower_fullDetector) <= 2 && (abs(phiOftower_fullDetector - cluster_phiOfTower_fullDetector) <= 2 or - abs(phiOftower_fullDetector - n_towers_Phi - cluster_phiOfTower_fullDetector) <= - 2)) { // The towers are within 3. Needs to stitch the two phi sides together + abs(phiOftower_fullDetector - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) { + // Remove the column outside of the L2 card: values (0,71), (23,26), (24,21), (47,50), (48,50), (71,2) if (!((cluster_phiOfTower_fullDetector == 0 && phiOftower_fullDetector == 71) or (cluster_phiOfTower_fullDetector == 23 && phiOftower_fullDetector == 26) or (cluster_phiOfTower_fullDetector == 24 && phiOftower_fullDetector == 21) or (cluster_phiOfTower_fullDetector == 47 && phiOftower_fullDetector == 50) or - (cluster_phiOfTower_fullDetector == n_links_GCTcard && phiOftower_fullDetector == 45) or - (cluster_phiOfTower_fullDetector == 71 && - phiOftower_fullDetector == 2))) { // Remove the column outside of the L2 card + (cluster_phiOfTower_fullDetector == 48 && phiOftower_fullDetector == 45) or + (cluster_phiOfTower_fullDetector == 71 && phiOftower_fullDetector == 2))) { isolation += ECAL_tower_L1Card[ll][mm][kk]; ntowers++; } } // Now do HCAL + // The towers are within 2. Needs to stitch the two phi sides together if (abs(etaOftower_fullDetector - cluster_etaOfTower_fullDetector) <= 2 && (abs(phiOftower_fullDetector - cluster_phiOfTower_fullDetector) <= 2 or - abs(phiOftower_fullDetector - n_towers_Phi - cluster_phiOfTower_fullDetector) <= - 2)) { // The towers are within 2. Needs to stitch the two phi sides together + abs(phiOftower_fullDetector - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) { hcal_nrj += HCAL_tower_L1Card[ll][mm][kk]; } } } } - cluster_list_L2[ii][jj].ciso = ((isolation) * (25.0 / ntowers)) / cluster_list_L2[ii][jj].cpt; - cluster_list_L2[ii][jj].chovere = hcal_nrj / cluster_list_L2[ii][jj].cpt; + cluster_list_L2[ii][jj].ciso_ = ((isolation) * (25.0 / ntowers)) / cluster_list_L2[ii][jj].cpt; + cluster_list_L2[ii][jj].chovere_ = hcal_nrj / cluster_list_L2[ii][jj].cpt; } } @@ -1041,29 +1035,29 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: for (int ii = 0; ii < n_towers_halfPhi; ++ii) { // The cluster list is still in the L1 like geometry for (unsigned int jj = 0; jj < unsigned(cluster_list_L2[ii].size()) && jj < n_clusters_4link; ++jj) { crystalID_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card] - [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ccrystalid; + [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ccrystalid_; towerID_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card] - [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ctowerid; + [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ctowerid_; energy_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card] [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cpt; brem_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card] - [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cbrem; + [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cbrem_; isolation_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card] - [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ciso; + [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ciso_; HE_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card] - [ii / n_clusters_4link] = cluster_list_L2[ii][jj].chovere; + [ii / n_clusters_4link] = cluster_list_L2[ii][jj].chovere_; showerShape_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card] - [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cshowershape; + [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cshowershape_; showerShapeLooseTk_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card] - [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cshowershapeloosetk; + [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cshowershapeloosetk_; photonShowerShape_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card] - [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cphotonshowershape; + [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cphotonshowershape_; } } - std::unique_ptr L1EGXtalClusters(new l1tp2::CaloCrystalClusterCollection); - std::unique_ptr > L1EGammas(new l1t::EGammaBxCollection); - std::unique_ptr L1CaloTowers(new l1tp2::CaloTowerCollection); + auto L1EGXtalClusters = std::make_unique(); + auto L1EGammas = std::make_unique(); + auto L1CaloTowers = std::make_unique(); // Fill the cluster collection and towers as well for (int ii = 0; ii < n_links_GCTcard; ++ii) { // 48 links @@ -1090,7 +1084,7 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: energy_cluster_L2Card[ii][jj][ll], HE_cluster_L2Card[ii][jj][ll], isolation_cluster_L2Card[ii][jj][ll], - centerhit.id, + centerhit.id(), -1000, float(brem_cluster_L2Card[ii][jj][ll]), -1000, @@ -1114,8 +1108,6 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: cluster.setExperimentalParams(params); L1EGXtalClusters->push_back(cluster); - // BXVector l1t::EGamma quality defined with respect to these WPs - // FIXME, need to defaul some of these to 0 I think... int standaloneWP = (int)(is_iso && is_ss); int looseL1TkMatchWP = (int)(is_looseTkiso && is_looseTkss); int photonWP = (int)(is_photon); @@ -1140,7 +1132,9 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: l1CaloTower.setTowerEta(getTowerEta_fromAbsoluteID(iEta_tower_L2Card[ii][jj][ll])); l1CaloTower.setTowerPhi(getTowerPhi_fromAbsoluteID(iPhi_tower_L2Card[ii][jj][ll])); // Some towers have incorrect eta/phi because that wasn't initialized in certain cases, fix these - if (l1CaloTower.towerEta() < -80. && l1CaloTower.towerPhi() < -90.) { + static float constexpr towerEtaUpperUnitialized = -80; + static float constexpr towerPhiUpperUnitialized = -90; + if (l1CaloTower.towerEta() < towerEtaUpperUnitialized && l1CaloTower.towerPhi() < towerPhiUpperUnitialized) { l1CaloTower.setTowerEta(l1t::CaloTools::towerEta(l1CaloTower.towerIEta())); l1CaloTower.setTowerPhi(l1t::CaloTools::towerPhi(l1CaloTower.towerIEta(), l1CaloTower.towerIPhi())); } @@ -1148,7 +1142,7 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: // Add L1EGs if they match in iEta / iPhi // L1EGs are already pT ordered, we will take the ID info for the leading one, but pT as the sum - for (auto l1eg : *L1EGXtalClusters) { + for (const auto l1eg : *L1EGXtalClusters) { if (l1eg.experimentalParam("TTiEta") != l1CaloTower.towerIEta()) continue; if (l1eg.experimentalParam("TTiPhi") != l1CaloTower.towerIPhi()) @@ -1157,8 +1151,9 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: int n_L1eg = l1CaloTower.nL1eg(); l1CaloTower.setNL1eg(n_L1eg++); l1CaloTower.setL1egTowerEt(l1CaloTower.l1egTowerEt() + l1eg.pt()); + // Don't record L1EG quality info for subleading L1EG if (l1CaloTower.nL1eg() > 1) - continue; // Don't record L1EG quality info for subleading L1EG + continue; l1CaloTower.setL1egTrkSS(l1eg.experimentalParam("trkMatchWP_showerShape")); l1CaloTower.setL1egTrkIso(l1eg.experimentalParam("trkMatchWP_isolation")); l1CaloTower.setL1egStandaloneSS(l1eg.experimentalParam("standaloneWP_showerShape")); @@ -1176,11 +1171,10 @@ void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm:: } bool L1EGCrystalClusterEmulatorProducer::passes_iso(float pt, float iso) { - if (pt < 80) { + if (pt < slideIsoPtThreshold) { if (!((a0_80 - a1_80 * pt) > iso)) return false; - } - if (pt >= 80) { + } else { if (iso > a0) return false; } diff --git a/L1Trigger/L1CaloTrigger/plugins/L1EGammaEEProducer.cc b/L1Trigger/L1CaloTrigger/plugins/L1EGammaEEProducer.cc index 943e5dc86dae0..85f0a1072db1e 100644 --- a/L1Trigger/L1CaloTrigger/plugins/L1EGammaEEProducer.cc +++ b/L1Trigger/L1CaloTrigger/plugins/L1EGammaEEProducer.cc @@ -9,17 +9,17 @@ #include "L1Trigger/L1CaloTrigger/interface/L1EGammaEECalibrator.h" #include "DataFormats/Math/interface/deltaPhi.h" -// we sort the clusters in pt -bool compare_cluster_pt(const l1t::HGCalMulticluster *cl1, const l1t::HGCalMulticluster *cl2); - -bool compare_cluster_pt(const l1t::HGCalMulticluster *cl1, const l1t::HGCalMulticluster *cl2) { - return cl1->pt() > cl2->pt(); -} +namespace l1tp2 { + // we sort the clusters in pt + bool compare_cluster_pt(const l1t::HGCalMulticluster *cl1, const l1t::HGCalMulticluster *cl2) { + return cl1->pt() > cl2->pt(); + } +}; // namespace l1tp2 int etaBin(const l1t::HGCalMulticluster *cl) { - float eta_min = 1.; - float eta_max = 4.; - unsigned n_eta_bins = 150; + static float constexpr eta_min = 1.; + static float constexpr eta_max = 4.; + static unsigned constexpr n_eta_bins = 150; int eta_bin = floor((std::abs(cl->eta()) - eta_min) / ((eta_max - eta_min) / n_eta_bins)); if (cl->eta() < 0) return -1 * eta_bin; // bin 0 doesn't exist @@ -27,9 +27,9 @@ int etaBin(const l1t::HGCalMulticluster *cl) { } int get_phi_bin(const l1t::HGCalMulticluster *cl) { - float phi_min = -M_PI; - float phi_max = M_PI; - unsigned n_phi_bins = 63; + static constexpr float phi_min = -M_PI; + static constexpr float phi_max = M_PI; + static constexpr unsigned n_phi_bins = 63; return floor(std::abs(reco::deltaPhi(cl->phi(), phi_min)) / ((phi_max - phi_min) / n_phi_bins)); } @@ -68,7 +68,6 @@ void L1EGammaEEProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSet // here we loop on the TPGs for (auto cl3d = multiclusters.begin(0); cl3d != multiclusters.end(0); cl3d++) { - // std::cout << "-- CL3D is EG: " << cl3d->hwQual() << " "<< cl3d->eta() <<" "<< cl3d->pt()<hwQual()) { if (cl3d->et() > minEt_) { int hw_quality = 1; // baseline EG ID passed @@ -107,13 +106,10 @@ void L1EGammaEEProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSet } } - // std::cout << "# of selected HGCAL multiclusters: " << selected_multiclusters.size() << std::endl; - - std::sort(selected_multiclusters.begin(), selected_multiclusters.end(), compare_cluster_pt); + std::sort(selected_multiclusters.begin(), selected_multiclusters.end(), l1tp2::compare_cluster_pt); std::set used_clusters; - for (auto cl3d : selected_multiclusters) { + for (const auto cl3d : selected_multiclusters) { if (used_clusters.find(cl3d) == used_clusters.end()) { - // reco::Candidate::LorentzVector mom = cl3d->p4(); float pt = cl3d->pt(); // we drop the Had component of the energy if (cl3d->hOverE() != -1) @@ -122,30 +118,19 @@ void L1EGammaEEProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSet reco::Candidate::PolarLorentzVector mom_eint( cl3d->iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM), cl3d->eta(), cl3d->phi(), 0.); - // cl3d->pt()/(1+cl3d->hOverE()) - // this is not yet used used_clusters.insert(cl3d); auto eta_phi_bin = get_eta_phi_bin(cl3d); - // std::cout << " cl pt: " << cl3d->pt() - // << " eta: " << cl3d->eta() - // << " phi: " << cl3d->phi() - // << " eta, phi bin: " << eta_phi_bin.first << "," << eta_phi_bin.second - // << std::endl; for (int eta_bin : {eta_phi_bin.first - 1, eta_phi_bin.first, eta_phi_bin.first + 1}) { for (int phi_bin : {eta_phi_bin.second - 1, eta_phi_bin.second, eta_phi_bin.second + 1}) { auto bucket = etaphi_bins.find(std::make_pair(eta_bin, phi_bin)); if (bucket != etaphi_bins.end()) { // this bucket is not empty - for (auto other_cl_ptr : bucket->second) { + for (const auto other_cl_ptr : bucket->second) { if (used_clusters.find(other_cl_ptr) == used_clusters.end()) { if (std::abs(other_cl_ptr->eta() - cl3d->eta()) < 0.02) { if (std::abs(reco::deltaPhi(other_cl_ptr->phi(), cl3d->phi())) < 0.1) { - // std::cout << "MERGE with cl pt: " << other_cl_ptr->pt() - // << " eta: " << other_cl_ptr->eta() - // << " phi: " << other_cl_ptr->phi() - // << std::endl; float pt_other = other_cl_ptr->pt(); if (other_cl_ptr->hOverE() != -1) pt_other = other_cl_ptr->pt() / (1 + other_cl_ptr->hOverE()); @@ -166,23 +151,19 @@ void L1EGammaEEProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSet float calib_factor = calibrator_.calibrationFactor(mom.pt(), mom.eta()); l1t::EGamma eg = l1t::EGamma(reco::Candidate::PolarLorentzVector(mom.pt() / calib_factor, mom.eta(), mom.phi(), 0.)); - // l1t::EGamma eg(mom); - // FIXME: full duplication with HWqual 2 eg.setHwQual(3); eg.setHwIso(1); l1EgammaBxCollection->push_back(0, eg); l1t::EGamma eg_emint_brec = l1t::EGamma(reco::Candidate::PolarLorentzVector(mom_eint.pt(), mom_eint.eta(), mom_eint.phi(), 0.)); - // l1t::EGamma eg(mom); - // FIXME: full duplication with HWqual 2 eg_emint_brec.setHwQual(5); eg_emint_brec.setHwIso(1); l1EgammaBxCollection->push_back(0, eg_emint_brec); } } - iEvent.put(std::move(l1EgammaBxCollection), "L1EGammaCollectionBXVWithCuts"); + iEvent.put(std::move(l1EgammaBxCollection)); } DEFINE_FWK_MODULE(L1EGammaEEProducer); diff --git a/L1Trigger/L1CaloTrigger/python/L1TowerCalibrationProducer_cfi.py b/L1Trigger/L1CaloTrigger/python/L1TowerCalibrationProducer_cfi.py index c65fb985d6638..dcadde6ad3a59 100644 --- a/L1Trigger/L1CaloTrigger/python/L1TowerCalibrationProducer_cfi.py +++ b/L1Trigger/L1CaloTrigger/python/L1TowerCalibrationProducer_cfi.py @@ -28,7 +28,6 @@ hfSF = cms.double(1.0), debug = cms.bool(False), skipCalibrations = cms.bool(False), - #debug = cms.bool(True), l1CaloTowers = cms.InputTag("L1EGammaClusterEmuProducer",), L1HgcalTowersInputTag = cms.InputTag("hgcalTowerProducer:HGCalTowerProcessor"), hcalDigis = cms.InputTag("simHcalTriggerPrimitiveDigis"), diff --git a/L1Trigger/L1CaloTrigger/python/l1EGammaEEProducer_cfi.py b/L1Trigger/L1CaloTrigger/python/l1EGammaEEProducer_cfi.py index 4bece082b7b88..161e8e8faa27e 100644 --- a/L1Trigger/L1CaloTrigger/python/l1EGammaEEProducer_cfi.py +++ b/L1Trigger/L1CaloTrigger/python/l1EGammaEEProducer_cfi.py @@ -1,5 +1,5 @@ import FWCore.ParameterSet.Config as cms l1EGammaEEProducer = cms.EDProducer("L1EGammaEEProducer", - calibrationConfig = cms.PSet(calirationFile = cms.FileInPath('L1Trigger/L1TCalorimeter/data/calib_ee_v1.json')), + calibrationConfig = cms.PSet(calibrationFile = cms.FileInPath('L1Trigger/L1TCalorimeter/data/calib_ee_v1.json')), Multiclusters=cms.InputTag('hgcalBackEndLayer2Producer:HGCalBackendLayer2Processor3DClustering')) diff --git a/L1Trigger/L1CaloTrigger/src/L1EGammaEECalibrator.cc b/L1Trigger/L1CaloTrigger/src/L1EGammaEECalibrator.cc index d8fb53c5113ad..ddd9ee590a12c 100644 --- a/L1Trigger/L1CaloTrigger/src/L1EGammaEECalibrator.cc +++ b/L1Trigger/L1CaloTrigger/src/L1EGammaEECalibrator.cc @@ -4,33 +4,33 @@ #include "boost/property_tree/ptree.hpp" #include "boost/property_tree/json_parser.hpp" #include -// #include -std::vector as_vector(boost::property_tree::ptree const& pt, boost::property_tree::ptree::key_type const& key) { - std::vector ret; - for (auto& item : pt.get_child(key)) - ret.push_back(item.second.get_value()); - return ret; -} +namespace l1tp2 { + std::vector as_vector(boost::property_tree::ptree const& pt, + boost::property_tree::ptree::key_type const& key) { + std::vector ret; + for (const auto& item : pt.get_child(key)) + ret.push_back(item.second.get_value()); + return ret; + } +}; // namespace l1tp2 L1EGammaEECalibrator::L1EGammaEECalibrator(const edm::ParameterSet& pset) { //read the JSON file and populate the eta - pt bins and the value container boost::property_tree::ptree calibration_map; - read_json(pset.getParameter("calirationFile").fullPath(), calibration_map); + read_json(pset.getParameter("calibrationFile").fullPath(), calibration_map); - auto eta_l = as_vector(calibration_map, "eta_l"); + auto eta_l = l1tp2::as_vector(calibration_map, "eta_l"); std::copy(eta_l.begin(), eta_l.end(), std::inserter(eta_bins, eta_bins.end())); - auto eta_h = as_vector(calibration_map, "eta_h"); + auto eta_h = l1tp2::as_vector(calibration_map, "eta_h"); eta_bins.insert(eta_h.back()); - // std::cout << "# of eta bins: " << eta_bins.size()-1 << std::endl; - auto pt_l = as_vector(calibration_map, "pt_l"); + auto pt_l = l1tp2::as_vector(calibration_map, "pt_l"); std::copy(pt_l.begin(), pt_l.end(), std::inserter(pt_bins, pt_bins.end())); - auto pt_h = as_vector(calibration_map, "pt_h"); + auto pt_h = l1tp2::as_vector(calibration_map, "pt_h"); pt_bins.insert(pt_h.back()); - // std::cout << "# of pt bins: " << pt_bins.size()-1 << std::endl; - auto calib_data = as_vector(calibration_map, "calib"); + auto calib_data = l1tp2::as_vector(calibration_map, "calib"); auto n_bins_eta = eta_bins.size(); auto n_bins_pt = pt_bins.size(); calib_factors.reserve(n_bins_eta * n_bins_pt); @@ -44,12 +44,10 @@ L1EGammaEECalibrator::L1EGammaEECalibrator(const edm::ParameterSet& pset) { int L1EGammaEECalibrator::bin(const std::set& container, float value) const { auto bin_l = container.upper_bound(value); - // std::cout << "value " << value << "lower boud " << *bin_l << " distance: " << std::distance(container.begin(), bin_l)-1 << std::endl; if (bin_l == container.end()) { // value not mapped to any bin return -1; } - // return bin_l - container.begin(); return std::distance(container.begin(), bin_l) - 1; } @@ -59,6 +57,5 @@ float L1EGammaEECalibrator::calibrationFactor(const float& pt, const float& eta) if (bin_eta == -1 || bin_pt == -1) return 1.; auto n_bins_pt = pt_bins.size(); - // std::cout << "pt: " << pt << " eta: " << eta << " ptbin: " << bin_pt << " etabin: " << bin_eta << " calib: " << calib_factors[(bin_eta*n_bins_pt)+bin_pt] << std::endl; return calib_factors[(bin_eta * n_bins_pt) + bin_pt]; } diff --git a/L1Trigger/L1CaloTrigger/src/ParametricCalibration.cc b/L1Trigger/L1CaloTrigger/src/ParametricCalibration.cc index 57eac0fb2f7d3..479bd2bcd4c21 100644 --- a/L1Trigger/L1CaloTrigger/src/ParametricCalibration.cc +++ b/L1Trigger/L1CaloTrigger/src/ParametricCalibration.cc @@ -1 +1,45 @@ #include "L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h" + +namespace l1tp2 { + ParametricCalibration::ParametricCalibration(const edm::ParameterSet& cpset) { + std::vector etaBins = cpset.getParameter>("etaBins"); + std::vector ptBins = cpset.getParameter>("ptBins"); + std::vector scale = cpset.getParameter>("scale"); + etas.insert(etas.end(), etaBins.begin(), etaBins.end()); + pts.insert(pts.end(), ptBins.begin(), ptBins.end()); + scales.insert(scales.end(), scale.begin(), scale.end()); + + std::vector ptMin = cpset.getParameter>("ptMin"); + if (ptMin.size() == etaBins.size()) { + ptMins.insert(ptMins.end(), ptMin.begin(), ptMin.end()); + } else if (ptMin.size() == 1) { + ptMins = std::vector(etaBins.size(), ptMin[0]); + } else { + throw cms::Exception("ParametricCalibration Configuration", "Ambiguous number of ptMin values in configuration"); + ; + } + + std::vector ptMax = cpset.getParameter>("ptMax"); + if (ptMax.size() == etaBins.size()) { + ptMaxs.insert(ptMaxs.end(), ptMax.begin(), ptMax.end()); + } else if (ptMax.size() == 1) { + ptMaxs = std::vector(etaBins.size(), ptMax[0]); + } else { + throw cms::Exception("ParametricCalibration Configuration", "Ambiguous number of ptMax values in configuration"); + ; + } + + if (pts.size() * etas.size() != scales.size()) + throw cms::Exception("Configuration", + "Bad number of calibration scales, pts.size() * etas.size() != scales.size()"); + } + // ------------ method fills 'descriptions' with the allowed parameters for the module ------------ + void ParametricCalibration::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.setComment(""); + desc.addUntracked>("ptMin", std::vector{}); + desc.addUntracked>("ptMax", std::vector{}); + descriptions.add("createIdealTkAlRecords", desc); + } + +}; // namespace l1tp2