File indexing completed on 2022-05-07 00:41:28
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021 #include <array>
0022 #include <iostream>
0023 #include <cmath>
0024
0025
0026 #include "FWCore/Framework/interface/stream/EDProducer.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/EventSetup.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030
0031 #include "CalibFormats/CaloTPG/interface/CaloTPGTranscoder.h"
0032 #include "CalibFormats/CaloTPG/interface/CaloTPGRecord.h"
0033 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0034 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0035 #include "Geometry/HcalTowerAlgo/interface/HcalTrigTowerGeometry.h"
0036 #include "FWCore/Framework/interface/MakerMacros.h"
0037 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0038 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0039
0040
0041 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0042
0043
0044 #include "DataFormats/HcalDigi/interface/HcalTriggerPrimitiveDigi.h"
0045
0046
0047 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloCrystalCluster.h"
0048 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloTower.h"
0049 #include "DataFormats/L1Trigger/interface/EGamma.h"
0050
0051 #include "L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h"
0052 #include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
0053 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0054
0055 static constexpr bool do_brem = true;
0056
0057 static constexpr int n_eta_bins = 2;
0058 static constexpr int n_borders_phi = 18;
0059 static constexpr int n_borders_eta = 18;
0060 static constexpr int n_clusters_max = 5;
0061 static constexpr int n_clusters_link = 3;
0062 static constexpr int n_clusters_4link = 4 * 3;
0063 static constexpr int n_crystals_towerEta = 5;
0064 static constexpr int n_crystals_towerPhi = 5;
0065 static constexpr int n_crystals_3towers = 3 * 5;
0066 static constexpr int n_towers_per_link = 17;
0067 static constexpr int n_clusters_per_link = 2;
0068 static constexpr int n_clusters_per_L1card = 8;
0069 static constexpr int n_towers_Eta = 34;
0070 static constexpr int n_towers_Phi = 72;
0071 static constexpr int n_towers_halfPhi = 36;
0072 static constexpr int n_links_card = 4;
0073 static constexpr int n_links_GCTcard = 48;
0074 static constexpr int n_GCTcards = 3;
0075 static constexpr float ECAL_eta_range = 1.4841;
0076 static constexpr float half_crystal_size = 0.00873;
0077 static constexpr float slideIsoPtThreshold = 80;
0078 static constexpr float plateau_ss = 130.0;
0079 static constexpr float a0_80 = 0.85, a1_80 = 0.0080, a0 = 0.21;
0080 static constexpr float b0 = 0.38, b1 = 1.9, b2 = 0.05;
0081 static constexpr float c0_ss = 0.94, c1_ss = 0.052, c2_ss = 0.044;
0082 static constexpr float d0 = 0.96, d1 = 0.0003;
0083 static constexpr float e0_looseTkss = 0.944, e1_looseTkss = 0.65, e2_looseTkss = 0.4;
0084 static constexpr float cut_500_MeV = 0.5;
0085
0086
0087
0088
0089 static constexpr int toweriEta_fromAbsoluteID_shift = 16;
0090
0091
0092
0093
0094 static constexpr int toweriPhi_fromAbsoluteID_shift_lowerHalf = 37;
0095 static constexpr int toweriPhi_fromAbsoluteID_shift_upperHalf = 35;
0096
0097 float getEta_fromL2LinkCardTowerCrystal(int link, int card, int tower, int crystal) {
0098 int etaID = n_crystals_towerEta * (n_towers_per_link * ((link / n_links_card) % 2) + (tower % n_towers_per_link)) +
0099 crystal % n_crystals_towerEta;
0100 float size_cell = 2 * ECAL_eta_range / (n_crystals_towerEta * n_towers_Eta);
0101 return etaID * size_cell - ECAL_eta_range + half_crystal_size;
0102 }
0103
0104 float getPhi_fromL2LinkCardTowerCrystal(int link, int card, int tower, int crystal) {
0105 int phiID = n_crystals_towerPhi * ((card * 24) + (n_links_card * (link / 8)) + (tower / n_towers_per_link)) +
0106 crystal / n_crystals_towerPhi;
0107 float size_cell = 2 * M_PI / (n_crystals_towerPhi * n_towers_Phi);
0108 return phiID * size_cell - M_PI + half_crystal_size;
0109 }
0110
0111 int getCrystal_etaID(float eta) {
0112 float size_cell = 2 * ECAL_eta_range / (n_crystals_towerEta * n_towers_Eta);
0113 int etaID = int((eta + ECAL_eta_range) / size_cell);
0114 return etaID;
0115 }
0116
0117 int convert_L2toL1_link(int link) { return link % n_links_card; }
0118
0119 int convert_L2toL1_tower(int tower) { return tower; }
0120
0121 int convert_L2toL1_card(int card, int link) { return card * n_clusters_4link + link / n_links_card; }
0122
0123 int getCrystal_phiID(float phi) {
0124 float size_cell = 2 * M_PI / (n_crystals_towerPhi * n_towers_Phi);
0125 int phiID = int((phi + M_PI) / size_cell);
0126 return phiID;
0127 }
0128
0129 int getTower_absoluteEtaID(float eta) {
0130 float size_cell = 2 * ECAL_eta_range / n_towers_Eta;
0131 int etaID = int((eta + ECAL_eta_range) / size_cell);
0132 return etaID;
0133 }
0134
0135 int getTower_absolutePhiID(float phi) {
0136 float size_cell = 2 * M_PI / n_towers_Phi;
0137 int phiID = int((phi + M_PI) / size_cell);
0138 return phiID;
0139 }
0140
0141 int getToweriEta_fromAbsoluteID(int id) {
0142 if (id < n_towers_per_link)
0143 return id - n_towers_per_link;
0144 else
0145 return id - toweriEta_fromAbsoluteID_shift;
0146 }
0147
0148 int getToweriPhi_fromAbsoluteID(int id) {
0149 if (id < n_towers_Phi / 2)
0150 return id + toweriPhi_fromAbsoluteID_shift_lowerHalf;
0151 else
0152 return id - toweriPhi_fromAbsoluteID_shift_upperHalf;
0153 }
0154
0155 float getTowerEta_fromAbsoluteID(int id) {
0156 float size_cell = 2 * ECAL_eta_range / n_towers_Eta;
0157 float eta = (id * size_cell) - ECAL_eta_range + 0.5 * size_cell;
0158 return eta;
0159 }
0160
0161 float getTowerPhi_fromAbsoluteID(int id) {
0162 float size_cell = 2 * M_PI / n_towers_Phi;
0163 float phi = (id * size_cell) - M_PI + 0.5 * size_cell;
0164 return phi;
0165 }
0166
0167 int getCrystalIDInTower(int etaID, int phiID) {
0168 return int(n_crystals_towerPhi * (phiID % n_crystals_towerPhi) + (etaID % n_crystals_towerEta));
0169 }
0170
0171 int get_towerEta_fromCardTowerInCard(int card, int towerincard) {
0172 return n_towers_per_link * (card % 2) + towerincard % n_towers_per_link;
0173 }
0174
0175 int get_towerPhi_fromCardTowerInCard(int card, int towerincard) {
0176 return 4 * (card / 2) + towerincard / n_towers_per_link;
0177 }
0178
0179 int get_towerEta_fromCardLinkTower(int card, int link, int tower) { return n_towers_per_link * (card % 2) + tower; }
0180
0181 int get_towerPhi_fromCardLinkTower(int card, int link, int tower) { return 4 * (card / 2) + link; }
0182
0183 int getTowerID(int etaID, int phiID) {
0184 return int(n_towers_per_link * ((phiID / n_crystals_towerPhi) % 4) +
0185 (etaID / n_crystals_towerEta) % n_towers_per_link);
0186 }
0187
0188 int getTower_phiID(int cluster_phiID) {
0189 return int((cluster_phiID / n_crystals_towerPhi) % 4);
0190 }
0191
0192 int getTower_etaID(int cluster_etaID) {
0193 return int((cluster_etaID / n_crystals_towerEta) % n_towers_per_link);
0194 }
0195
0196 int getEtaMax_card(int card) {
0197 int etamax = 0;
0198 if (card % 2 == 0)
0199 etamax = n_towers_per_link * n_crystals_towerEta - 1;
0200 else
0201 etamax = n_towers_Eta * n_crystals_towerEta - 1;
0202 return etamax;
0203 }
0204
0205 int getEtaMin_card(int card) {
0206 int etamin = 0;
0207 if (card % 2 == 0)
0208 etamin = 0 * n_crystals_towerEta;
0209 else
0210 etamin = n_towers_per_link * n_crystals_towerEta;
0211 return etamin;
0212 }
0213
0214 int getPhiMax_card(int card) {
0215 int phimax = ((card / 2) + 1) * 4 * n_crystals_towerPhi - 1;
0216 return phimax;
0217 }
0218
0219 int getPhiMin_card(int card) {
0220 int phimin = (card / 2) * 4 * n_crystals_towerPhi;
0221 return phimin;
0222 }
0223
0224 class L1EGCrystalClusterEmulatorProducer : public edm::stream::EDProducer<> {
0225 public:
0226 explicit L1EGCrystalClusterEmulatorProducer(const edm::ParameterSet&);
0227 ~L1EGCrystalClusterEmulatorProducer() override;
0228
0229 private:
0230 void produce(edm::Event&, const edm::EventSetup&) override;
0231 bool passes_ss(float pt, float ss);
0232 bool passes_photon(float pt, float pss);
0233 bool passes_iso(float pt, float iso);
0234 bool passes_looseTkss(float pt, float ss);
0235 bool passes_looseTkiso(float pt, float iso);
0236 float get_calibrate(float uncorr);
0237
0238 edm::EDGetTokenT<EcalEBTrigPrimDigiCollection> ecalTPEBToken_;
0239 edm::EDGetTokenT<edm::SortedCollection<HcalTriggerPrimitiveDigi> > hcalTPToken_;
0240 edm::ESGetToken<CaloTPGTranscoder, CaloTPGRecord> decoderTag_;
0241
0242 l1tp2::ParametricCalibration calib_;
0243
0244 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryTag_;
0245 const CaloSubdetectorGeometry* ebGeometry;
0246 const CaloSubdetectorGeometry* hbGeometry;
0247 edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> hbTopologyTag_;
0248 const HcalTopology* hcTopology_;
0249
0250 struct mycluster {
0251 float c2x2_;
0252 float c2x5_;
0253 float c5x5_;
0254 int cshowershape_;
0255 int cshowershapeloosetk_;
0256 float cvalueshowershape_;
0257 int cphotonshowershape_;
0258 float cpt;
0259 int cbrem_;
0260 float cWeightedEta_;
0261 float cWeightedPhi_;
0262 float ciso_;
0263 float chovere_;
0264 float craweta_;
0265 float crawphi_;
0266 float chcal_;
0267 float ceta_;
0268 float cphi_;
0269 int ccrystalid_;
0270 int cinsidecrystalid_;
0271 int ctowerid_;
0272 };
0273
0274 class SimpleCaloHit {
0275 private:
0276 float pt_ = 0;
0277 float energy_ = 0.;
0278 bool isEndcapHit_ = false;
0279 bool stale_ = false;
0280 bool used_ = false;
0281 GlobalVector position_;
0282 HcalDetId id_hcal_;
0283 EBDetId id_;
0284
0285 public:
0286
0287 inline void setPt() { pt_ = (position_.mag2() > 0) ? energy_ * sin(position_.theta()) : 0; };
0288 inline void setEnergy(float et) { energy_ = et / sin(position_.theta()); };
0289 inline void setIsEndcapHit(bool isEC) { isEndcapHit_ = isEC; };
0290 inline void setUsed(bool isUsed) { used_ = isUsed; };
0291 inline void setPosition(const GlobalVector& pos) { position_ = pos; };
0292 inline void setIdHcal(const HcalDetId& idhcal) { id_hcal_ = idhcal; };
0293 inline void setId(const EBDetId& id) { id_ = id; };
0294
0295 inline float pt() const { return pt_; };
0296 inline float energy() const { return energy_; };
0297 inline bool isEndcapHit() const { return isEndcapHit_; };
0298 inline bool used() const { return used_; };
0299 inline const GlobalVector& position() const { return position_; };
0300 inline const EBDetId& id() const { return id_; };
0301
0302 inline float deta(SimpleCaloHit& other) const { return position_.eta() - other.position().eta(); };
0303 int dieta(SimpleCaloHit& other) const {
0304 if (isEndcapHit_ || other.isEndcapHit())
0305 return 9999;
0306 if (id_.ieta() * other.id().ieta() > 0)
0307 return id_.ieta() - other.id().ieta();
0308 return id_.ieta() - other.id().ieta() - 1;
0309 };
0310 inline float dphi(SimpleCaloHit& other) const {
0311 return reco::deltaPhi(static_cast<float>(position_.phi()), static_cast<float>(other.position().phi()));
0312 };
0313 int diphi(SimpleCaloHit& other) const {
0314 if (isEndcapHit_ || other.isEndcapHit())
0315 return 9999;
0316
0317 static constexpr int PI = 180;
0318 int result = id().iphi() - other.id().iphi();
0319 while (result > PI)
0320 result -= 2 * PI;
0321 while (result <= -PI)
0322 result += 2 * PI;
0323 return result;
0324 };
0325 float distanceTo(SimpleCaloHit& other) const {
0326
0327
0328 return (position() - other.position()).mag();
0329 };
0330 bool operator==(SimpleCaloHit& other) const {
0331 return (id_ == other.id() && position() == other.position() && energy_ == other.energy() &&
0332 isEndcapHit_ == other.isEndcapHit());
0333 };
0334 };
0335 };
0336
0337 L1EGCrystalClusterEmulatorProducer::L1EGCrystalClusterEmulatorProducer(const edm::ParameterSet& iConfig)
0338 : ecalTPEBToken_(consumes<EcalEBTrigPrimDigiCollection>(iConfig.getParameter<edm::InputTag>("ecalTPEB"))),
0339 hcalTPToken_(
0340 consumes<edm::SortedCollection<HcalTriggerPrimitiveDigi> >(iConfig.getParameter<edm::InputTag>("hcalTP"))),
0341 decoderTag_(esConsumes<CaloTPGTranscoder, CaloTPGRecord>(edm::ESInputTag("", ""))),
0342 calib_(iConfig.getParameter<edm::ParameterSet>("calib")),
0343 caloGeometryTag_(esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag("", ""))),
0344 hbTopologyTag_(esConsumes<HcalTopology, HcalRecNumberingRecord>(edm::ESInputTag("", ""))) {
0345 produces<l1tp2::CaloCrystalClusterCollection>();
0346 produces<BXVector<l1t::EGamma> >();
0347 produces<l1tp2::CaloTowerCollection>("L1CaloTowerCollection");
0348 }
0349
0350 L1EGCrystalClusterEmulatorProducer::~L1EGCrystalClusterEmulatorProducer() {}
0351
0352 void L1EGCrystalClusterEmulatorProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0353 using namespace edm;
0354
0355 edm::Handle<EcalEBTrigPrimDigiCollection> pcalohits;
0356 iEvent.getByToken(ecalTPEBToken_, pcalohits);
0357
0358 const auto& caloGeometry = iSetup.getData(caloGeometryTag_);
0359 ebGeometry = caloGeometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0360 hbGeometry = caloGeometry.getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
0361 const auto& hbTopology = iSetup.getData(hbTopologyTag_);
0362 hcTopology_ = &hbTopology;
0363 HcalTrigTowerGeometry theTrigTowerGeometry(hcTopology_);
0364
0365 const auto& decoder = iSetup.getData(decoderTag_);
0366
0367
0368
0369
0370
0371
0372 iEvent.getByToken(ecalTPEBToken_, pcalohits);
0373 std::vector<SimpleCaloHit> ecalhits;
0374
0375 for (const auto& hit : *pcalohits.product()) {
0376 if (hit.encodedEt() > 0)
0377 {
0378
0379 float et = hit.encodedEt() / 8.;
0380 if (et < cut_500_MeV)
0381 continue;
0382
0383 auto cell = ebGeometry->getGeometry(hit.id());
0384
0385 SimpleCaloHit ehit;
0386 ehit.setId(hit.id());
0387 ehit.setPosition(GlobalVector(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z()));
0388 ehit.setEnergy(et);
0389 ehit.setPt();
0390 ecalhits.push_back(ehit);
0391 }
0392 }
0393
0394
0395 std::vector<SimpleCaloHit> hcalhits;
0396 edm::Handle<edm::SortedCollection<HcalTriggerPrimitiveDigi> > hbhecoll;
0397 iEvent.getByToken(hcalTPToken_, hbhecoll);
0398 for (const auto& hit : *hbhecoll.product()) {
0399 float et = decoder.hcaletValue(hit.id(), hit.t0());
0400 if (et <= 0)
0401 continue;
0402 if (!(hcTopology_->validHT(hit.id()))) {
0403 LogError("L1EGCrystalClusterEmulatorProducer")
0404 << " -- Hcal hit DetID not present in HCAL Geom: " << hit.id() << std::endl;
0405 throw cms::Exception("L1EGCrystalClusterEmulatorProducer");
0406 continue;
0407 }
0408 const std::vector<HcalDetId>& hcId = theTrigTowerGeometry.detIds(hit.id());
0409 if (hcId.empty()) {
0410 LogError("L1EGCrystalClusterEmulatorProducer")
0411 << "Cannot find any HCalDetId corresponding to " << hit.id() << std::endl;
0412 throw cms::Exception("L1EGCrystalClusterEmulatorProducer");
0413 continue;
0414 }
0415 if (hcId[0].subdetId() > 1)
0416 continue;
0417 GlobalVector hcal_tp_position = GlobalVector(0., 0., 0.);
0418 for (const auto& hcId_i : hcId) {
0419 if (hcId_i.subdetId() > 1)
0420 continue;
0421 auto cell = hbGeometry->getGeometry(hcId_i);
0422 if (cell == nullptr)
0423 continue;
0424 GlobalVector tmpVector = GlobalVector(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z());
0425 hcal_tp_position = tmpVector;
0426 break;
0427 }
0428 SimpleCaloHit hhit;
0429 hhit.setId(hit.id());
0430 hhit.setIdHcal(hit.id());
0431 hhit.setPosition(hcal_tp_position);
0432 hhit.setEnergy(et);
0433 hhit.setPt();
0434 hcalhits.push_back(hhit);
0435 }
0436
0437
0438
0439
0440
0441
0442
0443 float ECAL_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
0444 float HCAL_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
0445 int iEta_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
0446 int iPhi_tower_L1Card[n_links_card][n_towers_per_link][n_towers_halfPhi];
0447
0448 float energy_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
0449
0450 int brem_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
0451 int towerID_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
0452 int crystalID_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
0453 int showerShape_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
0454 int showerShapeLooseTk_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
0455 int photonShowerShape_cluster_L1Card[n_links_card][n_clusters_link][n_towers_halfPhi];
0456
0457 for (int ii = 0; ii < n_links_card; ++ii) {
0458 for (int jj = 0; jj < n_towers_per_link; ++jj) {
0459 for (int ll = 0; ll < n_towers_halfPhi; ++ll) {
0460 ECAL_tower_L1Card[ii][jj][ll] = 0;
0461 HCAL_tower_L1Card[ii][jj][ll] = 0;
0462 iPhi_tower_L1Card[ii][jj][ll] = -999;
0463 iEta_tower_L1Card[ii][jj][ll] = -999;
0464 }
0465 }
0466 }
0467 for (int ii = 0; ii < n_links_card; ++ii) {
0468 for (int jj = 0; jj < n_clusters_link; ++jj) {
0469 for (int ll = 0; ll < n_towers_halfPhi; ++ll) {
0470 energy_cluster_L1Card[ii][jj][ll] = 0;
0471 brem_cluster_L1Card[ii][jj][ll] = 0;
0472 towerID_cluster_L1Card[ii][jj][ll] = 0;
0473 crystalID_cluster_L1Card[ii][jj][ll] = 0;
0474 }
0475 }
0476 }
0477
0478
0479 vector<mycluster> cluster_list[n_towers_halfPhi];
0480
0481 vector<mycluster> cluster_list_merged[n_towers_halfPhi];
0482
0483 for (int cc = 0; cc < n_towers_halfPhi; ++cc) {
0484
0485 for (int nregion = 0; nregion <= n_clusters_max; ++nregion) {
0486 int nclusters = 0;
0487 bool build_cluster = true;
0488
0489
0490 while (nclusters < n_clusters_max && build_cluster) {
0491 build_cluster = false;
0492 SimpleCaloHit centerhit;
0493
0494 for (const auto& hit : ecalhits) {
0495 if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
0496 getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
0497 getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
0498 getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) &&
0499
0500 getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) &&
0501 getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion &&
0502 !hit.used() && hit.pt() >= 1.0 && hit.pt() > centerhit.pt())
0503 {
0504
0505 centerhit = hit;
0506 build_cluster = true;
0507 }
0508 }
0509 if (build_cluster)
0510 nclusters++;
0511
0512
0513 if (build_cluster && nclusters > 0 && nclusters <= n_clusters_max) {
0514 mycluster mc1;
0515 mc1.cpt = 0.0;
0516 mc1.cWeightedEta_ = 0.0;
0517 mc1.cWeightedPhi_ = 0.0;
0518 float leftlobe = 0;
0519 float rightlobe = 0;
0520 float e5x5 = 0;
0521 float n5x5 = 0;
0522 float e2x5_1 = 0;
0523 float n2x5_1 = 0;
0524 float e2x5_2 = 0;
0525 float n2x5_2 = 0;
0526 float e2x2_1 = 0;
0527 float n2x2_1 = 0;
0528 float e2x2_2 = 0;
0529 float n2x2_2 = 0;
0530 float e2x2_3 = 0;
0531 float n2x2_3 = 0;
0532 float e2x2_4 = 0;
0533 float n2x2_4 = 0;
0534 for (auto& hit : ecalhits) {
0535 if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
0536 getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
0537 getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
0538 getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && hit.pt() > 0 &&
0539 getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) &&
0540 getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion) {
0541 if (abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) > 2 && hit.diphi(centerhit) <= 7) {
0542 rightlobe += hit.pt();
0543 }
0544 if (abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) < -2 && hit.diphi(centerhit) >= -7) {
0545 leftlobe += hit.pt();
0546 }
0547 if (abs(hit.dieta(centerhit)) <= 2 && abs(hit.diphi(centerhit)) <= 2) {
0548 e5x5 += hit.energy();
0549 n5x5++;
0550 }
0551 if ((hit.dieta(centerhit) == 1 or hit.dieta(centerhit) == 0) &&
0552 (hit.diphi(centerhit) == 1 or hit.diphi(centerhit) == 0)) {
0553 e2x2_1 += hit.energy();
0554 n2x2_1++;
0555 }
0556 if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) &&
0557 (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == 1)) {
0558 e2x2_2 += hit.energy();
0559 n2x2_2++;
0560 }
0561 if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == 1) &&
0562 (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == -1)) {
0563 e2x2_3 += hit.energy();
0564 n2x2_3++;
0565 }
0566 if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) &&
0567 (hit.diphi(centerhit) == 0 or hit.diphi(centerhit) == -1)) {
0568 e2x2_4 += hit.energy();
0569 n2x2_4++;
0570 }
0571 if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == 1) && abs(hit.diphi(centerhit)) <= 2) {
0572 e2x5_1 += hit.energy();
0573 n2x5_1++;
0574 }
0575 if ((hit.dieta(centerhit) == 0 or hit.dieta(centerhit) == -1) && abs(hit.diphi(centerhit)) <= 2) {
0576 e2x5_2 += hit.energy();
0577 n2x5_2++;
0578 }
0579 }
0580 if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
0581 getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
0582 getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
0583 getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && !hit.used() && hit.pt() > 0 &&
0584 abs(hit.dieta(centerhit)) <= 1 && abs(hit.diphi(centerhit)) <= 2 &&
0585 getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) &&
0586 getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion) {
0587
0588 hit.setUsed(true);
0589 mc1.cpt += hit.pt();
0590 mc1.cWeightedEta_ += float(hit.pt()) * float(hit.position().eta());
0591 mc1.cWeightedPhi_ = mc1.cWeightedPhi_ + (float(hit.pt()) * float(hit.position().phi()));
0592 }
0593 }
0594 if (do_brem && (rightlobe > 0.10 * mc1.cpt or leftlobe > 0.10 * mc1.cpt)) {
0595 for (auto& hit : ecalhits) {
0596 if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
0597 getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
0598 getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
0599 getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && hit.pt() > 0 &&
0600 getCrystal_etaID(hit.position().eta()) < getEtaMin_card(cc) + n_crystals_3towers * (nregion + 1) &&
0601 getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) + n_crystals_3towers * nregion &&
0602 !hit.used()) {
0603 if (rightlobe > 0.10 * mc1.cpt && (leftlobe < 0.10 * mc1.cpt or rightlobe > leftlobe) &&
0604 abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) > 2 && hit.diphi(centerhit) <= 7) {
0605 mc1.cpt += hit.pt();
0606 hit.setUsed(true);
0607 mc1.cbrem_ = 1;
0608 }
0609 if (leftlobe > 0.10 * mc1.cpt && (rightlobe < 0.10 * mc1.cpt or leftlobe >= rightlobe) &&
0610 abs(hit.dieta(centerhit)) <= 1 && hit.diphi(centerhit) < -2 && hit.diphi(centerhit) >= -7) {
0611 mc1.cpt += hit.pt();
0612 hit.setUsed(true);
0613 mc1.cbrem_ = 1;
0614 }
0615 }
0616 }
0617 }
0618 mc1.c5x5_ = e5x5;
0619 mc1.c2x5_ = max(e2x5_1, e2x5_2);
0620 mc1.c2x2_ = e2x2_1;
0621 if (e2x2_2 > mc1.c2x2_)
0622 mc1.c2x2_ = e2x2_2;
0623 if (e2x2_3 > mc1.c2x2_)
0624 mc1.c2x2_ = e2x2_3;
0625 if (e2x2_4 > mc1.c2x2_)
0626 mc1.c2x2_ = e2x2_4;
0627 mc1.cWeightedEta_ = mc1.cWeightedEta_ / mc1.cpt;
0628 mc1.cWeightedPhi_ = mc1.cWeightedPhi_ / mc1.cpt;
0629 mc1.ceta_ = getCrystal_etaID(centerhit.position().eta());
0630 mc1.cphi_ = getCrystal_phiID(centerhit.position().phi());
0631 mc1.crawphi_ = centerhit.position().phi();
0632 mc1.craweta_ = centerhit.position().eta();
0633 cluster_list[cc].push_back(mc1);
0634 }
0635 }
0636 }
0637 std::sort(begin(cluster_list[cc]), end(cluster_list[cc]), [](mycluster a, mycluster b) { return a.cpt > b.cpt; });
0638
0639
0640 for (unsigned int jj = 0; jj < unsigned(cluster_list[cc].size()); ++jj) {
0641 for (unsigned int kk = jj + 1; kk < unsigned(cluster_list[cc].size()); ++kk) {
0642 if (std::abs(cluster_list[cc][jj].ceta_ - cluster_list[cc][kk].ceta_) < 2 &&
0643 std::abs(cluster_list[cc][jj].cphi_ - cluster_list[cc][kk].cphi_) < 2) {
0644 if (cluster_list[cc][kk].cpt > cluster_list[cc][jj].cpt) {
0645 cluster_list[cc][kk].cpt += cluster_list[cc][jj].cpt;
0646 cluster_list[cc][kk].c5x5_ += cluster_list[cc][jj].c5x5_;
0647 cluster_list[cc][kk].c2x5_ += cluster_list[cc][jj].c2x5_;
0648 cluster_list[cc][jj].cpt = 0;
0649 cluster_list[cc][jj].c5x5_ = 0;
0650 cluster_list[cc][jj].c2x5_ = 0;
0651 cluster_list[cc][jj].c2x2_ = 0;
0652 } else {
0653 cluster_list[cc][jj].cpt += cluster_list[cc][kk].cpt;
0654 cluster_list[cc][jj].c5x5_ += cluster_list[cc][kk].c5x5_;
0655 cluster_list[cc][jj].c2x5_ += cluster_list[cc][kk].c2x5_;
0656 cluster_list[cc][kk].cpt = 0;
0657 cluster_list[cc][kk].c2x2_ = 0;
0658 cluster_list[cc][kk].c2x5_ = 0;
0659 cluster_list[cc][kk].c5x5_ = 0;
0660 }
0661 }
0662 }
0663 if (cluster_list[cc][jj].cpt > 0) {
0664 cluster_list[cc][jj].cpt =
0665 cluster_list[cc][jj].cpt *
0666 calib_(cluster_list[cc][jj].cpt,
0667 std::abs(cluster_list[cc][jj].craweta_));
0668 cluster_list_merged[cc].push_back(cluster_list[cc][jj]);
0669 }
0670 }
0671 std::sort(begin(cluster_list_merged[cc]), end(cluster_list_merged[cc]), [](mycluster a, mycluster b) {
0672 return a.cpt > b.cpt;
0673 });
0674
0675
0676 for (unsigned int jj = 0; jj < unsigned(cluster_list_merged[cc].size()) && jj < n_clusters_4link; ++jj) {
0677 crystalID_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] =
0678 getCrystalIDInTower(cluster_list_merged[cc][jj].ceta_, cluster_list_merged[cc][jj].cphi_);
0679 towerID_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] =
0680 getTowerID(cluster_list_merged[cc][jj].ceta_, cluster_list_merged[cc][jj].cphi_);
0681 energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = cluster_list_merged[cc][jj].cpt;
0682 brem_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = cluster_list_merged[cc][jj].cbrem_;
0683 if (passes_ss(cluster_list_merged[cc][jj].cpt,
0684 cluster_list_merged[cc][jj].c2x5_ / cluster_list_merged[cc][jj].c5x5_))
0685 showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1;
0686 else
0687 showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0;
0688 if (passes_looseTkss(cluster_list_merged[cc][jj].cpt,
0689 cluster_list_merged[cc][jj].c2x5_ / cluster_list_merged[cc][jj].c5x5_))
0690 showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1;
0691 else
0692 showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0;
0693 if (passes_photon(cluster_list_merged[cc][jj].cpt,
0694 cluster_list_merged[cc][jj].c2x2_ / cluster_list_merged[cc][jj].c2x5_))
0695 photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 1;
0696 else
0697 photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][cc] = 0;
0698 }
0699
0700
0701 for (const auto& hit : ecalhits) {
0702 if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
0703 getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
0704 getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
0705 getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) &&
0706 !hit.used()) {
0707 for (int jj = 0; jj < n_links_card; ++jj) {
0708 if ((getCrystal_phiID(hit.position().phi()) / n_crystals_towerPhi) % 4 == jj) {
0709 for (int ii = 0; ii < n_towers_per_link; ++ii) {
0710
0711 if ((getCrystal_etaID(hit.position().eta()) / n_crystals_towerEta) % n_towers_per_link == ii) {
0712 ECAL_tower_L1Card[jj][ii][cc] += hit.pt() * calib_(0, std::abs(hit.position().eta()));
0713 iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(hit.position().eta());
0714 iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(hit.position().phi());
0715 }
0716 }
0717 }
0718 }
0719
0720 static constexpr float tower_width = 0.0873;
0721 for (int jj = 0; jj < n_links_card; ++jj) {
0722 for (int ii = 0; ii < n_towers_per_link; ++ii) {
0723 float phi = getPhiMin_card(cc) * tower_width / n_crystals_towerPhi - M_PI + (jj + 0.5) * tower_width;
0724 float eta = getEtaMin_card(cc) * tower_width / n_crystals_towerEta - n_towers_per_link * tower_width +
0725 (ii + 0.5) * tower_width;
0726 iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(eta);
0727 iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(phi);
0728 }
0729 }
0730
0731 }
0732 }
0733
0734
0735 for (const auto& hit : hcalhits) {
0736 if (getCrystal_phiID(hit.position().phi()) <= getPhiMax_card(cc) &&
0737 getCrystal_phiID(hit.position().phi()) >= getPhiMin_card(cc) &&
0738 getCrystal_etaID(hit.position().eta()) <= getEtaMax_card(cc) &&
0739 getCrystal_etaID(hit.position().eta()) >= getEtaMin_card(cc) && hit.pt() > 0) {
0740 for (int jj = 0; jj < n_links_card; ++jj) {
0741 if ((getCrystal_phiID(hit.position().phi()) / n_crystals_towerPhi) % n_links_card == jj) {
0742 for (int ii = 0; ii < n_towers_per_link; ++ii) {
0743 if ((getCrystal_etaID(hit.position().eta()) / n_crystals_towerEta) % n_towers_per_link == ii) {
0744 HCAL_tower_L1Card[jj][ii][cc] += hit.pt();
0745 iEta_tower_L1Card[jj][ii][cc] = getTower_absoluteEtaID(hit.position().eta());
0746 iPhi_tower_L1Card[jj][ii][cc] = getTower_absolutePhiID(hit.position().phi());
0747 }
0748 }
0749 }
0750 }
0751 }
0752 }
0753
0754
0755 for (unsigned int kk = n_clusters_4link; kk < cluster_list_merged[cc].size(); ++kk) {
0756 if (cluster_list_merged[cc][kk].cpt > 0) {
0757 ECAL_tower_L1Card[getTower_phiID(cluster_list_merged[cc][kk].cphi_)]
0758 [getTower_etaID(cluster_list_merged[cc][kk].ceta_)][cc] += cluster_list_merged[cc][kk].cpt;
0759 }
0760 }
0761 }
0762
0763
0764
0765
0766
0767
0768 float HCAL_tower_L2Card[n_links_GCTcard][n_towers_per_link]
0769 [n_GCTcards];
0770 float ECAL_tower_L2Card[n_links_GCTcard][n_towers_per_link][n_GCTcards];
0771 int iEta_tower_L2Card[n_links_GCTcard][n_towers_per_link][n_GCTcards];
0772 int iPhi_tower_L2Card[n_links_GCTcard][n_towers_per_link][n_GCTcards];
0773 float energy_cluster_L2Card[n_links_GCTcard][n_clusters_per_link]
0774 [n_GCTcards];
0775 float brem_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0776 int towerID_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0777 int crystalID_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0778 float isolation_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0779 float HE_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0780 int showerShape_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0781 int showerShapeLooseTk_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0782 int photonShowerShape_cluster_L2Card[n_links_GCTcard][n_clusters_per_link][n_GCTcards];
0783
0784 for (int ii = 0; ii < n_links_GCTcard; ++ii) {
0785 for (int jj = 0; jj < n_towers_per_link; ++jj) {
0786 for (int ll = 0; ll < n_GCTcards; ++ll) {
0787 HCAL_tower_L2Card[ii][jj][ll] = 0;
0788 ECAL_tower_L2Card[ii][jj][ll] = 0;
0789 iEta_tower_L2Card[ii][jj][ll] = 0;
0790 iPhi_tower_L2Card[ii][jj][ll] = 0;
0791 }
0792 }
0793 }
0794 for (int ii = 0; ii < n_links_GCTcard; ++ii) {
0795 for (int jj = 0; jj < n_clusters_per_link; ++jj) {
0796 for (int ll = 0; ll < n_GCTcards; ++ll) {
0797 energy_cluster_L2Card[ii][jj][ll] = 0;
0798 brem_cluster_L2Card[ii][jj][ll] = 0;
0799 towerID_cluster_L2Card[ii][jj][ll] = 0;
0800 crystalID_cluster_L2Card[ii][jj][ll] = 0;
0801 isolation_cluster_L2Card[ii][jj][ll] = 0;
0802 HE_cluster_L2Card[ii][jj][ll] = 0;
0803 photonShowerShape_cluster_L2Card[ii][jj][ll] = 0;
0804 showerShape_cluster_L2Card[ii][jj][ll] = 0;
0805 showerShapeLooseTk_cluster_L2Card[ii][jj][ll] = 0;
0806 }
0807 }
0808 }
0809
0810
0811 vector<mycluster> cluster_list_L2[n_towers_halfPhi];
0812
0813
0814 for (int ii = 0; ii < n_borders_phi; ++ii) {
0815 for (int jj = 0; jj < n_eta_bins; ++jj) {
0816 int card_left = 2 * ii + jj;
0817 int card_right = 2 * ii + jj + 2;
0818 if (card_right > 35)
0819 card_right = card_right - n_towers_halfPhi;
0820 for (int kk = 0; kk < n_clusters_4link; ++kk) {
0821 if (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 50 &&
0822 crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 19 &&
0823 energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 0) {
0824 for (int ll = 0; ll < n_clusters_4link; ++ll) {
0825 if (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < n_towers_per_link &&
0826 crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < 5 &&
0827 std::abs(
0828 5 * (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right]) % n_towers_per_link +
0829 crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] % 5 -
0830 5 * (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left]) % n_towers_per_link -
0831 crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] % 5) < 2) {
0832 if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] >
0833 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right]) {
0834 energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] +=
0835 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right];
0836 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] = 0;
0837 }
0838 else {
0839 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] +=
0840 energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left];
0841 energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] = 0;
0842 }
0843 }
0844 }
0845 }
0846 }
0847 }
0848 }
0849
0850
0851 if (do_brem) {
0852 for (int ii = 0; ii < n_borders_phi; ++ii) {
0853 for (int jj = 0; jj < n_eta_bins; ++jj) {
0854 int card_left = 2 * ii + jj;
0855 int card_right = 2 * ii + jj + 2;
0856 if (card_right > 35)
0857 card_right = card_right - n_towers_halfPhi;
0858
0859 for (int kk = 0; kk < n_clusters_4link; ++kk) {
0860
0861 if (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 50 &&
0862 energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] > 0) {
0863 for (int ll = 0; ll < n_clusters_4link; ++ll) {
0864
0865 if (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < n_towers_per_link &&
0866 std::abs(5 * (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right]) %
0867 n_towers_per_link +
0868 crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] % 5 -
0869 5 * (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left]) %
0870 n_towers_per_link -
0871 crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] % 5) <= 1) {
0872
0873 if (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] < n_towers_per_link &&
0874 std::abs(5 + crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] / 5 -
0875 (crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] / 5)) <= 5) {
0876 if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] >
0877 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] &&
0878 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] >
0879 0.10 * energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left]) {
0880 energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] +=
0881 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right];
0882 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] = 0;
0883 brem_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] = 1;
0884 }
0885
0886 else if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_right] >=
0887 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_left] &&
0888 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_left] >
0889 0.10 * energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_right]) {
0890 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] +=
0891 energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left];
0892 energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_left] = 0;
0893 brem_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_right] = 1;
0894 }
0895 }
0896 }
0897 }
0898 }
0899 }
0900 }
0901 }
0902 }
0903
0904
0905 for (int ii = 0; ii < n_borders_eta; ++ii) {
0906 int card_bottom = 2 * ii;
0907 int card_top = 2 * ii + 1;
0908 for (int kk = 0; kk < n_clusters_4link; ++kk) {
0909 if (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] % n_towers_per_link == 16 &&
0910 crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] % 5 == n_links_card &&
0911 energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] >
0912 0) {
0913 for (int ll = 0; ll < n_clusters_4link; ++ll) {
0914 if (std::abs(
0915 5 * (towerID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] / n_towers_per_link) +
0916 crystalID_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] / 5 -
0917 5 * (towerID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] / n_towers_per_link) -
0918 crystalID_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] / 5) < 2) {
0919 if (energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] >
0920 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_bottom]) {
0921 energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] +=
0922 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top];
0923 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] = 0;
0924 } else {
0925 energy_cluster_L1Card[ll % n_links_card][ll / n_links_card][card_top] +=
0926 energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom];
0927 energy_cluster_L1Card[kk % n_links_card][kk / n_links_card][card_bottom] = 0;
0928 }
0929 }
0930 }
0931 }
0932 }
0933 }
0934
0935
0936 for (int ii = 0; ii < n_towers_halfPhi; ++ii) {
0937 for (int jj = 0; jj < n_clusters_4link; ++jj) {
0938 if (energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii] > 0) {
0939 mycluster mc1;
0940 mc1.cpt = energy_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
0941 mc1.cbrem_ = brem_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
0942 mc1.ctowerid_ = towerID_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
0943 mc1.ccrystalid_ = crystalID_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
0944 mc1.cshowershape_ = showerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
0945 mc1.cshowershapeloosetk_ = showerShapeLooseTk_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
0946 mc1.cphotonshowershape_ = photonShowerShape_cluster_L1Card[jj % n_links_card][jj / n_links_card][ii];
0947 cluster_list_L2[ii].push_back(mc1);
0948 }
0949 }
0950 std::sort(
0951 begin(cluster_list_L2[ii]), end(cluster_list_L2[ii]), [](mycluster a, mycluster b) { return a.cpt > b.cpt; });
0952 }
0953
0954
0955 for (int ii = 0; ii < n_towers_halfPhi; ++ii) {
0956 for (unsigned int jj = n_clusters_per_L1card; jj < n_clusters_4link && jj < cluster_list_L2[ii].size(); ++jj) {
0957 if (cluster_list_L2[ii][jj].cpt > 0) {
0958 ECAL_tower_L1Card[cluster_list_L2[ii][jj].ctowerid_ / n_towers_per_link]
0959 [cluster_list_L2[ii][jj].ctowerid_ % n_towers_per_link][ii] += cluster_list_L2[ii][jj].cpt;
0960 cluster_list_L2[ii][jj].cpt = 0;
0961 cluster_list_L2[ii][jj].ctowerid_ = 0;
0962 cluster_list_L2[ii][jj].ccrystalid_ = 0;
0963 }
0964 }
0965 }
0966
0967
0968 for (int ii = 0; ii < n_towers_halfPhi; ++ii) {
0969 for (unsigned int jj = 0; jj < n_clusters_per_L1card && jj < cluster_list_L2[ii].size(); ++jj) {
0970 int cluster_etaOfTower_fullDetector = get_towerEta_fromCardTowerInCard(ii, cluster_list_L2[ii][jj].ctowerid_);
0971 int cluster_phiOfTower_fullDetector = get_towerPhi_fromCardTowerInCard(ii, cluster_list_L2[ii][jj].ctowerid_);
0972 float hcal_nrj = 0.0;
0973 float isolation = 0.0;
0974 int ntowers = 0;
0975 for (int iii = 0; iii < n_towers_halfPhi; ++iii) {
0976 for (unsigned int jjj = 0; jjj < n_clusters_per_L1card && jjj < cluster_list_L2[iii].size(); ++jjj) {
0977 if (!(iii == ii && jjj == jj)) {
0978 int cluster2_eta = get_towerEta_fromCardTowerInCard(iii, cluster_list_L2[iii][jjj].ctowerid_);
0979 int cluster2_phi = get_towerPhi_fromCardTowerInCard(iii, cluster_list_L2[iii][jjj].ctowerid_);
0980 if (abs(cluster2_eta - cluster_etaOfTower_fullDetector) <= 2 &&
0981 (abs(cluster2_phi - cluster_phiOfTower_fullDetector) <= 2 or
0982 abs(cluster2_phi - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) {
0983 isolation += cluster_list_L2[iii][jjj].cpt;
0984 }
0985 }
0986 }
0987 }
0988 for (int kk = 0; kk < n_towers_halfPhi; ++kk) {
0989 for (int ll = 0; ll < n_links_card; ++ll) {
0990 for (int mm = 0; mm < n_towers_per_link; ++mm) {
0991 int etaOftower_fullDetector = get_towerEta_fromCardLinkTower(kk, ll, mm);
0992 int phiOftower_fullDetector = get_towerPhi_fromCardLinkTower(kk, ll, mm);
0993
0994
0995 if (abs(etaOftower_fullDetector - cluster_etaOfTower_fullDetector) <= 2 &&
0996 (abs(phiOftower_fullDetector - cluster_phiOfTower_fullDetector) <= 2 or
0997 abs(phiOftower_fullDetector - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) {
0998
0999 if (!((cluster_phiOfTower_fullDetector == 0 && phiOftower_fullDetector == 71) or
1000 (cluster_phiOfTower_fullDetector == 23 && phiOftower_fullDetector == 26) or
1001 (cluster_phiOfTower_fullDetector == 24 && phiOftower_fullDetector == 21) or
1002 (cluster_phiOfTower_fullDetector == 47 && phiOftower_fullDetector == 50) or
1003 (cluster_phiOfTower_fullDetector == 48 && phiOftower_fullDetector == 45) or
1004 (cluster_phiOfTower_fullDetector == 71 && phiOftower_fullDetector == 2))) {
1005 isolation += ECAL_tower_L1Card[ll][mm][kk];
1006 ntowers++;
1007 }
1008 }
1009
1010
1011 if (abs(etaOftower_fullDetector - cluster_etaOfTower_fullDetector) <= 2 &&
1012 (abs(phiOftower_fullDetector - cluster_phiOfTower_fullDetector) <= 2 or
1013 abs(phiOftower_fullDetector - n_towers_Phi - cluster_phiOfTower_fullDetector) <= 2)) {
1014 hcal_nrj += HCAL_tower_L1Card[ll][mm][kk];
1015 }
1016 }
1017 }
1018 }
1019 cluster_list_L2[ii][jj].ciso_ = ((isolation) * (25.0 / ntowers)) / cluster_list_L2[ii][jj].cpt;
1020 cluster_list_L2[ii][jj].chovere_ = hcal_nrj / cluster_list_L2[ii][jj].cpt;
1021 }
1022 }
1023
1024
1025
1026 for (int ii = 0; ii < n_links_GCTcard; ++ii) {
1027 for (int jj = 0; jj < n_towers_per_link; ++jj) {
1028 for (int ll = 0; ll < 3; ++ll) {
1029 ECAL_tower_L2Card[ii][jj][ll] =
1030 ECAL_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
1031 HCAL_tower_L2Card[ii][jj][ll] =
1032 HCAL_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
1033 iEta_tower_L2Card[ii][jj][ll] =
1034 iEta_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
1035 iPhi_tower_L2Card[ii][jj][ll] =
1036 iPhi_tower_L1Card[convert_L2toL1_link(ii)][convert_L2toL1_tower(jj)][convert_L2toL1_card(ll, ii)];
1037 }
1038 }
1039 }
1040
1041
1042 for (int ii = 0; ii < n_towers_halfPhi; ++ii) {
1043 for (unsigned int jj = 0; jj < unsigned(cluster_list_L2[ii].size()) && jj < n_clusters_per_L1card; ++jj) {
1044 crystalID_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1045 [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ccrystalid_;
1046 towerID_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1047 [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ctowerid_;
1048 energy_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1049 [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cpt;
1050 brem_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1051 [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cbrem_;
1052 isolation_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1053 [ii / n_clusters_4link] = cluster_list_L2[ii][jj].ciso_;
1054 HE_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1055 [ii / n_clusters_4link] = cluster_list_L2[ii][jj].chovere_;
1056 showerShape_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1057 [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cshowershape_;
1058 showerShapeLooseTk_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1059 [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cshowershapeloosetk_;
1060 photonShowerShape_cluster_L2Card[n_links_card * (ii % n_clusters_4link) + jj % n_links_card][jj / n_links_card]
1061 [ii / n_clusters_4link] = cluster_list_L2[ii][jj].cphotonshowershape_;
1062 }
1063 }
1064
1065 auto L1EGXtalClusters = std::make_unique<l1tp2::CaloCrystalClusterCollection>();
1066 auto L1EGammas = std::make_unique<l1t::EGammaBxCollection>();
1067 auto L1CaloTowers = std::make_unique<l1tp2::CaloTowerCollection>();
1068
1069
1070 for (int ii = 0; ii < n_links_GCTcard; ++ii) {
1071 for (int ll = 0; ll < n_GCTcards; ++ll) {
1072
1073 for (int jj = 0; jj < 2; ++jj) {
1074 if (energy_cluster_L2Card[ii][jj][ll] > 0.45) {
1075 reco::Candidate::PolarLorentzVector p4calibrated(
1076 energy_cluster_L2Card[ii][jj][ll],
1077 getEta_fromL2LinkCardTowerCrystal(
1078 ii, ll, towerID_cluster_L2Card[ii][jj][ll], crystalID_cluster_L2Card[ii][jj][ll]),
1079 getPhi_fromL2LinkCardTowerCrystal(
1080 ii, ll, towerID_cluster_L2Card[ii][jj][ll], crystalID_cluster_L2Card[ii][jj][ll]),
1081 0.);
1082 SimpleCaloHit centerhit;
1083 bool is_iso = passes_iso(energy_cluster_L2Card[ii][jj][ll], isolation_cluster_L2Card[ii][jj][ll]);
1084 bool is_looseTkiso =
1085 passes_looseTkiso(energy_cluster_L2Card[ii][jj][ll], isolation_cluster_L2Card[ii][jj][ll]);
1086 bool is_ss = (showerShape_cluster_L2Card[ii][jj][ll] == 1);
1087 bool is_photon = (photonShowerShape_cluster_L2Card[ii][jj][ll] == 1) && is_ss && is_iso;
1088 bool is_looseTkss = (showerShapeLooseTk_cluster_L2Card[ii][jj][ll] == 1);
1089
1090 l1tp2::CaloCrystalCluster cluster(p4calibrated,
1091 energy_cluster_L2Card[ii][jj][ll],
1092 HE_cluster_L2Card[ii][jj][ll],
1093 isolation_cluster_L2Card[ii][jj][ll],
1094 centerhit.id(),
1095 -1000,
1096 float(brem_cluster_L2Card[ii][jj][ll]),
1097 -1000,
1098 -1000,
1099 energy_cluster_L2Card[ii][jj][ll],
1100 -1,
1101 is_iso&& is_ss,
1102 is_iso&& is_ss,
1103 is_photon,
1104 is_iso&& is_ss,
1105 is_looseTkiso&& is_looseTkss,
1106 is_iso&& is_ss);
1107
1108 std::map<std::string, float> params;
1109 params["standaloneWP_showerShape"] = is_ss;
1110 params["standaloneWP_isolation"] = is_iso;
1111 params["trkMatchWP_showerShape"] = is_looseTkss;
1112 params["trkMatchWP_isolation"] = is_looseTkiso;
1113 params["TTiEta"] = getToweriEta_fromAbsoluteID(getTower_absoluteEtaID(p4calibrated.eta()));
1114 params["TTiPhi"] = getToweriPhi_fromAbsoluteID(getTower_absolutePhiID(p4calibrated.phi()));
1115 cluster.setExperimentalParams(params);
1116 L1EGXtalClusters->push_back(cluster);
1117
1118 int standaloneWP = (int)(is_iso && is_ss);
1119 int looseL1TkMatchWP = (int)(is_looseTkiso && is_looseTkss);
1120 int photonWP = (int)(is_photon);
1121 int quality =
1122 (standaloneWP * std::pow(2, 0)) + (looseL1TkMatchWP * std::pow(2, 1)) + (photonWP * std::pow(2, 2));
1123 L1EGammas->push_back(
1124 0, l1t::EGamma(p4calibrated, p4calibrated.pt(), p4calibrated.eta(), p4calibrated.phi(), quality, 1));
1125 }
1126 }
1127 }
1128 }
1129
1130 for (int ii = 0; ii < n_links_GCTcard; ++ii) {
1131 for (int ll = 0; ll < n_GCTcards; ++ll) {
1132
1133 for (int jj = 0; jj < n_towers_per_link; ++jj) {
1134 l1tp2::CaloTower l1CaloTower;
1135 l1CaloTower.setEcalTowerEt(ECAL_tower_L2Card[ii][jj][ll]);
1136 l1CaloTower.setHcalTowerEt(HCAL_tower_L2Card[ii][jj][ll]);
1137 l1CaloTower.setTowerIEta(getToweriEta_fromAbsoluteID(iEta_tower_L2Card[ii][jj][ll]));
1138 l1CaloTower.setTowerIPhi(getToweriPhi_fromAbsoluteID(iPhi_tower_L2Card[ii][jj][ll]));
1139 l1CaloTower.setTowerEta(getTowerEta_fromAbsoluteID(iEta_tower_L2Card[ii][jj][ll]));
1140 l1CaloTower.setTowerPhi(getTowerPhi_fromAbsoluteID(iPhi_tower_L2Card[ii][jj][ll]));
1141
1142 static float constexpr towerEtaUpperUnitialized = -80;
1143 static float constexpr towerPhiUpperUnitialized = -90;
1144 if (l1CaloTower.towerEta() < towerEtaUpperUnitialized && l1CaloTower.towerPhi() < towerPhiUpperUnitialized) {
1145 l1CaloTower.setTowerEta(l1t::CaloTools::towerEta(l1CaloTower.towerIEta()));
1146 l1CaloTower.setTowerPhi(l1t::CaloTools::towerPhi(l1CaloTower.towerIEta(), l1CaloTower.towerIPhi()));
1147 }
1148 l1CaloTower.setIsBarrel(true);
1149
1150
1151
1152 for (const auto& l1eg : *L1EGXtalClusters) {
1153 if (l1eg.experimentalParam("TTiEta") != l1CaloTower.towerIEta())
1154 continue;
1155 if (l1eg.experimentalParam("TTiPhi") != l1CaloTower.towerIPhi())
1156 continue;
1157
1158 int n_L1eg = l1CaloTower.nL1eg();
1159 l1CaloTower.setNL1eg(n_L1eg++);
1160 l1CaloTower.setL1egTowerEt(l1CaloTower.l1egTowerEt() + l1eg.pt());
1161
1162 if (l1CaloTower.nL1eg() > 1)
1163 continue;
1164 l1CaloTower.setL1egTrkSS(l1eg.experimentalParam("trkMatchWP_showerShape"));
1165 l1CaloTower.setL1egTrkIso(l1eg.experimentalParam("trkMatchWP_isolation"));
1166 l1CaloTower.setL1egStandaloneSS(l1eg.experimentalParam("standaloneWP_showerShape"));
1167 l1CaloTower.setL1egStandaloneIso(l1eg.experimentalParam("standaloneWP_isolation"));
1168 }
1169
1170 L1CaloTowers->push_back(l1CaloTower);
1171 }
1172 }
1173 }
1174
1175 iEvent.put(std::move(L1EGXtalClusters));
1176 iEvent.put(std::move(L1EGammas));
1177 iEvent.put(std::move(L1CaloTowers), "L1CaloTowerCollection");
1178 }
1179
1180 bool L1EGCrystalClusterEmulatorProducer::passes_iso(float pt, float iso) {
1181 bool is_iso = true;
1182 if (pt < slideIsoPtThreshold) {
1183 if (!((a0_80 - a1_80 * pt) > iso))
1184 is_iso = false;
1185 } else {
1186 if (iso > a0)
1187 is_iso = false;
1188 }
1189 if (pt > plateau_ss)
1190 is_iso = true;
1191 return is_iso;
1192 }
1193
1194 bool L1EGCrystalClusterEmulatorProducer::passes_looseTkiso(float pt, float iso) {
1195 bool is_iso = (b0 + b1 * std::exp(-b2 * pt) > iso);
1196 if (pt > plateau_ss)
1197 is_iso = true;
1198 return is_iso;
1199 }
1200
1201 bool L1EGCrystalClusterEmulatorProducer::passes_ss(float pt, float ss) {
1202 bool is_ss = ((c0_ss + c1_ss * std::exp(-c2_ss * pt)) <= ss);
1203 if (pt > plateau_ss)
1204 is_ss = true;
1205 return is_ss;
1206 }
1207
1208 bool L1EGCrystalClusterEmulatorProducer::passes_photon(float pt, float pss) {
1209 bool is_ss = (pss > d0 - d1 * pt);
1210 if (pt > plateau_ss)
1211 is_ss = true;
1212 return is_ss;
1213 }
1214
1215 bool L1EGCrystalClusterEmulatorProducer::passes_looseTkss(float pt, float ss) {
1216 bool is_ss = ((e0_looseTkss - e1_looseTkss * std::exp(-e2_looseTkss * pt)) <= ss);
1217 if (pt > plateau_ss)
1218 is_ss = true;
1219 return is_ss;
1220 }
1221
1222
1223 DEFINE_FWK_MODULE(L1EGCrystalClusterEmulatorProducer);