L1CaloTrigger/plugins/Phase2L1CaloEGammaEmulator.cc

0001 /*
0002  * Description: Phase 2 RCT and GCT emulator
0003  */
0004
0005 // system include files
0006 #include <ap_int.h>
0007 #include <array>
0008 #include <cmath>
0009 // #include <cstdint>
0010 #include <cstdlib>  // for rand
0011 #include <iostream>
0012 #include <fstream>
0013 #include <memory>
0014
0015 // user include files
0016 #include "FWCore/Framework/interface/stream/EDProducer.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/ESHandle.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0023
0024 #include "CalibFormats/CaloTPG/interface/CaloTPGTranscoder.h"
0025 #include "CalibFormats/CaloTPG/interface/CaloTPGRecord.h"
0026 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0027 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0028 #include "Geometry/HcalTowerAlgo/interface/HcalTrigTowerGeometry.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0031 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0032
0033 // ECAL TPs
0034 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0035
0036 // HCAL TPs
0037 #include "DataFormats/HcalDigi/interface/HcalTriggerPrimitiveDigi.h"
0038
0039 // Output tower collection
0040 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloCrystalCluster.h"
0041 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloTower.h"
0042 #include "DataFormats/L1TCalorimeterPhase2/interface/DigitizedClusterCorrelator.h"
0043 #include "DataFormats/L1TCalorimeterPhase2/interface/DigitizedTowerCorrelator.h"
0044 #include "DataFormats/L1TCalorimeterPhase2/interface/DigitizedClusterGT.h"
0045
0046 #include "DataFormats/L1Trigger/interface/BXVector.h"
0047 #include "DataFormats/L1Trigger/interface/EGamma.h"
0048
0049 #include "L1Trigger/L1CaloTrigger/interface/ParametricCalibration.h"
0050 #include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
0051
0052 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0053
0054 #include "L1Trigger/L1CaloTrigger/interface/Phase2L1CaloEGammaUtils.h"
0055 #include "L1Trigger/L1CaloTrigger/interface/Phase2L1RCT.h"
0056 #include "L1Trigger/L1CaloTrigger/interface/Phase2L1GCT.h"
0057
0058 // Declare the Phase2L1CaloEGammaEmulator class and its methods
0059
0060 class Phase2L1CaloEGammaEmulator : public edm::stream::EDProducer<> {
0061 public:
0062   explicit Phase2L1CaloEGammaEmulator(const edm::ParameterSet&);
0063   ~Phase2L1CaloEGammaEmulator() override = default;
0064
0065   static void fillDescriptions(edm::ConfigurationDescriptions&);
0066
0067 private:
0068   void produce(edm::Event&, const edm::EventSetup&) override;
0069
0070   edm::EDGetTokenT<EcalEBTrigPrimDigiCollection> ecalTPEBToken_;
0071   edm::EDGetTokenT<edm::SortedCollection<HcalTriggerPrimitiveDigi>> hcalTPToken_;
0072   edm::ESGetToken<CaloTPGTranscoder, CaloTPGRecord> decoderTag_;
0073
0074   l1tp2::ParametricCalibration calib_;
0075
0076   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryTag_;
0077   const CaloSubdetectorGeometry* ebGeometry;
0078   const CaloSubdetectorGeometry* hbGeometry;
0079   edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> hbTopologyTag_;
0080   const HcalTopology* hcTopology_;
0081 };
0082
0083 //////////////////////////////////////////////////////////////////////////
0084
0085 // Phase2L1CaloEGammaEmulator initializer, destructor, and produce methods
0086
0087 Phase2L1CaloEGammaEmulator::Phase2L1CaloEGammaEmulator(const edm::ParameterSet& iConfig)
0088     : ecalTPEBToken_(consumes<EcalEBTrigPrimDigiCollection>(iConfig.getParameter<edm::InputTag>("ecalTPEB"))),
0089       hcalTPToken_(
0090           consumes<edm::SortedCollection<HcalTriggerPrimitiveDigi>>(iConfig.getParameter<edm::InputTag>("hcalTP"))),
0091       decoderTag_(esConsumes<CaloTPGTranscoder, CaloTPGRecord>(edm::ESInputTag("", ""))),
0092       calib_(iConfig.getParameter<edm::ParameterSet>("calib")),
0093       caloGeometryTag_(esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag("", ""))),
0094       hbTopologyTag_(esConsumes<HcalTopology, HcalRecNumberingRecord>(edm::ESInputTag("", ""))) {
0095   produces<l1tp2::CaloCrystalClusterCollection>("RCTClusters");
0096   produces<l1tp2::CaloCrystalClusterCollection>("GCTClusters");
0097   produces<l1tp2::CaloTowerCollection>("RCTTowers");
0098   produces<l1tp2::CaloTowerCollection>("GCTTowers");
0099   produces<l1tp2::CaloTowerCollection>("GCTFullTowers");
0100   produces<BXVector<l1t::EGamma>>("GCTEGammas");
0101   produces<l1tp2::DigitizedClusterCorrelatorCollection>("GCTDigitizedClusterToCorrelator");
0102   produces<l1tp2::DigitizedTowerCorrelatorCollection>("GCTDigitizedTowerToCorrelator");
0103   produces<l1tp2::DigitizedClusterGTCollection>("GCTDigitizedClusterToGT");
0104 }
0105
0106 void Phase2L1CaloEGammaEmulator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0107   using namespace edm;
0108
0109   // Detector geometry
0110   const auto& caloGeometry = iSetup.getData(caloGeometryTag_);
0111   ebGeometry = caloGeometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0112   hbGeometry = caloGeometry.getSubdetectorGeometry(DetId::Hcal, HcalBarrel);
0113   const auto& hbTopology = iSetup.getData(hbTopologyTag_);
0114   hcTopology_ = &hbTopology;
0115   HcalTrigTowerGeometry theTrigTowerGeometry(hcTopology_);
0116
0117   const auto& decoder = iSetup.getData(decoderTag_);
0118
0119   //***************************************************//
0120   // Declare RCT output collections
0121   //***************************************************//
0122
0123   auto L1EGXtalClusters = std::make_unique<l1tp2::CaloCrystalClusterCollection>();
0124   auto L1CaloTowers = std::make_unique<l1tp2::CaloTowerCollection>();
0125
0126   //***************************************************//
0127   // Get the ECAL hits
0128   //***************************************************//
0129   edm::Handle<EcalEBTrigPrimDigiCollection> pcalohits;
0130   iEvent.getByToken(ecalTPEBToken_, pcalohits);
0131
0132   std::vector<p2eg::SimpleCaloHit> ecalhits;
0133
0134   for (const auto& hit : *pcalohits.product()) {
0135     if (hit.encodedEt() > 0)  // hit.encodedEt() returns an int corresponding to 2x the crystal Et
0136     {
0137       // Et is 10 bit, by keeping the ADC saturation Et at 120 GeV it means that you have to multiply by 0.125 (input LSB)
0138       float et = hit.encodedEt() * 0.125;
0139       if (et < p2eg::cut_500_MeV) {
0140         continue;  // Reject hits with < 500 MeV ET
0141       }
0142
0143       // Get cell coordinates and info
0144       auto cell = ebGeometry->getGeometry(hit.id());
0145
0146       p2eg::SimpleCaloHit ehit;
0147       ehit.setId(hit.id());
0148       ehit.setPosition(GlobalVector(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z()));
0149       ehit.setEnergy(et);
0150       ehit.setEt_uint(
0151           (ap_uint<10>)hit.encodedEt() >>
0152           2);  // also save the uint Et, this is to convert between 0.125 (in MC production) and 0.5 (in firmware based code)
0153       ehit.setPt();
0154       ecalhits.push_back(ehit);
0155     }
0156   }
0157
0158   //***************************************************//
0159   // Get the HCAL hits
0160   //***************************************************//
0161   std::vector<p2eg::SimpleCaloHit> hcalhits;
0162   edm::Handle<edm::SortedCollection<HcalTriggerPrimitiveDigi>> hbhecoll;
0163   iEvent.getByToken(hcalTPToken_, hbhecoll);
0164
0165   for (const auto& hit : *hbhecoll.product()) {
0166     float et = decoder.hcaletValue(hit.id(), hit.t0());
0167     ap_uint<10> encodedEt = hit.t0().compressedEt();
0168     // same thing as SOI_compressedEt() in HcalTriggerPrimitiveDigi.h///
0169     if (et <= 0)
0170       continue;
0171
0172     if (!(hcTopology_->validHT(hit.id()))) {
0173       LogError("Phase2L1CaloEGammaEmulator")
0174           << " -- Hcal hit DetID not present in HCAL Geom: " << hit.id() << std::endl;
0175       throw cms::Exception("Phase2L1CaloEGammaEmulator");
0176       continue;
0177     }
0178     const std::vector<HcalDetId>& hcId = theTrigTowerGeometry.detIds(hit.id());
0179     if (hcId.empty()) {
0180       LogError("Phase2L1CaloEGammaEmulator") << "Cannot find any HCalDetId corresponding to " << hit.id() << std::endl;
0181       throw cms::Exception("Phase2L1CaloEGammaEmulator");
0182       continue;
0183     }
0184     if (hcId[0].subdetId() > 1) {
0185       continue;
0186     }
0187     GlobalVector hcal_tp_position = GlobalVector(0., 0., 0.);
0188     for (const auto& hcId_i : hcId) {
0189       if (hcId_i.subdetId() > 1) {
0190         continue;
0191       }
0192       // get the first HCAL TP/ cell
0193       auto cell = hbGeometry->getGeometry(hcId_i);
0194       if (cell == nullptr) {
0195         continue;
0196       }
0197       GlobalVector tmpVector = GlobalVector(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z());
0198       hcal_tp_position = tmpVector;
0199
0200       break;
0201     }
0202     p2eg::SimpleCaloHit hhit;
0203     hhit.setId(hit.id());
0204     hhit.setIdHcal(hit.id());
0205     hhit.setPosition(hcal_tp_position);
0206     hhit.setEnergy(et);
0207     hhit.setPt();
0208     hhit.setEt_uint(encodedEt);
0209     hcalhits.push_back(hhit);
0210   }
0211
0212   //***************************************************//
0213   // Initialize necessary arrays for tower and clusters
0214   //***************************************************//
0215
0216   // L1 Outputs definition: Arrays that use firmware convention for indexing
0217   p2eg::tower_t towerHCALCard
0218       [p2eg::n_towers_cardEta][p2eg::n_towers_cardPhi]
0219       [p2eg::n_towers_halfPhi];  // 17x4x36 array (not to be confused with the 12x1 array of ap_uints, towerEtHCAL
0220   p2eg::tower_t towerECALCard[p2eg::n_towers_cardEta][p2eg::n_towers_cardPhi][p2eg::n_towers_halfPhi];
0221   // There is one vector of clusters per card (up to 12 clusters before stitching across ECAL regions)
0222   std::vector<p2eg::Cluster> cluster_list[p2eg::n_towers_halfPhi];
0223   // After merging/stitching the clusters, we only take the 8 highest pt per card
0224   std::vector<p2eg::Cluster> cluster_list_merged[p2eg::n_towers_halfPhi];
0225
0226   //***************************************************//
0227   // Fill RCT ECAL regions with ECAL hits
0228   //***************************************************//
0229   for (int cc = 0; cc < p2eg::n_towers_halfPhi; ++cc) {  // Loop over 36 L1 cards
0230
0231     p2eg::card rctCard;
0232     rctCard.setIdx(cc);
0233
0234     for (const auto& hit : ecalhits) {
0235       // Check if the hit is in cards 0-35
0236       if (hit.isInCard(cc)) {
0237         // Get the crystal eta and phi, relative to the bottom left corner of the card
0238         // (0 up to 17*5, 0 up to 4*5)
0239         int local_iEta = hit.crystalLocaliEta(cc);
0240         int local_iPhi = hit.crystalLocaliPhi(cc);
0241
0242         // Region number (0-5) depends only on the crystal iEta in the card
0243         int regionNumber = p2eg::getRegionNumber(local_iEta);
0244
0245         // Tower eta and phi index inside the card (17x4)
0246         int inCard_tower_iEta = local_iEta / p2eg::CRYSTALS_IN_TOWER_ETA;
0247         int inCard_tower_iPhi = local_iPhi / p2eg::CRYSTALS_IN_TOWER_PHI;
0248
0249         // Tower eta and phi index inside the region (3x4)
0250         int inRegion_tower_iEta = inCard_tower_iEta % p2eg::TOWER_IN_ETA;
0251         int inRegion_tower_iPhi = inCard_tower_iPhi % p2eg::TOWER_IN_PHI;
0252
0253         // Crystal eta and phi index inside the 3x4 region (15x20)
0254         int inRegion_crystal_iEta = local_iEta % (p2eg::TOWER_IN_ETA * p2eg::CRYSTALS_IN_TOWER_ETA);
0255         int inRegion_crystal_iPhi = local_iPhi;
0256
0257         // Crystal eta and phi index inside the tower (5x5)
0258         int inLink_crystal_iEta = (inRegion_crystal_iEta % p2eg::CRYSTALS_IN_TOWER_ETA);
0259         int inLink_crystal_iPhi = (inRegion_crystal_iPhi % p2eg::CRYSTALS_IN_TOWER_PHI);
0260
0261         // Add the crystal energy to the rctCard
0262         p2eg::region3x4& myRegion = rctCard.getRegion3x4(regionNumber);
0263         p2eg::linkECAL& myLink = myRegion.getLinkECAL(inRegion_tower_iEta, inRegion_tower_iPhi);
0264         myLink.addCrystalE(inLink_crystal_iEta, inLink_crystal_iPhi, hit.et_uint());
0265       }
0266     }
0267
0268     //***************************************************//
0269     // Build RCT towers from HCAL hits
0270     //***************************************************//
0271     for (const auto& hit : hcalhits) {
0272       if (hit.isInCard(cc) && hit.pt() > 0) {
0273         // Get crystal eta and phi, relative to the bottom left corner of the card
0274         // (0 up to 17*5, 0 up to 4*5)
0275         int local_iEta = hit.crystalLocaliEta(cc);
0276         int local_iPhi = hit.crystalLocaliPhi(cc);
0277
0278         // Region (0-5) the hit falls into
0279         int regionNumber = p2eg::getRegionNumber(local_iEta);
0280
0281         // Tower eta and phi index inside the card (17x4)
0282         int inCard_tower_iEta = int(local_iEta / p2eg::CRYSTALS_IN_TOWER_ETA);
0283         int inCard_tower_iPhi = int(local_iPhi / p2eg::CRYSTALS_IN_TOWER_PHI);
0284
0285         // Tower eta and phi index inside the region (3x4)
0286         int inRegion_tower_iEta = inCard_tower_iEta % p2eg::TOWER_IN_ETA;
0287         int inRegion_tower_iPhi = inCard_tower_iPhi % p2eg::TOWER_IN_PHI;
0288
0289         // Access the right HCAL region and tower and increment the ET
0290         p2eg::towers3x4& myTowers3x4 = rctCard.getTowers3x4(regionNumber);
0291         p2eg::towerHCAL& myTower = myTowers3x4.getTowerHCAL(inRegion_tower_iEta, inRegion_tower_iPhi);
0292         myTower.addEt(hit.et_uint());
0293       }
0294     }
0295
0296     //***************************************************//
0297     // Make clusters in each ECAL region independently
0298     //***************************************************//
0299     for (int idxRegion = 0; idxRegion < p2eg::N_REGIONS_PER_CARD; idxRegion++) {
0300       // ECAL crystals array
0301       p2eg::crystal temporary[p2eg::CRYSTAL_IN_ETA][p2eg::CRYSTAL_IN_PHI];
0302       // HCAL towers in 3x4 region
0303       ap_uint<12> towerEtHCAL[p2eg::TOWER_IN_ETA * p2eg::TOWER_IN_PHI];
0304
0305       p2eg::region3x4& myRegion = rctCard.getRegion3x4(idxRegion);
0306       p2eg::towers3x4& myTowers = rctCard.getTowers3x4(idxRegion);
0307
0308       // In each 3x4 region, loop through the links (one link = one tower)
0309       for (int iLinkEta = 0; iLinkEta < p2eg::TOWER_IN_ETA; iLinkEta++) {
0310         for (int iLinkPhi = 0; iLinkPhi < p2eg::TOWER_IN_PHI; iLinkPhi++) {
0311           // Get the ECAL link (one link per tower)
0312           p2eg::linkECAL& myLink = myRegion.getLinkECAL(iLinkEta, iLinkPhi);
0313
0314           // We have an array of 3x4 links/towers, each link/tower is 5x5 in crystals. We need to convert this to a 15x20 of crystals
0315           int ref_iEta = (iLinkEta * p2eg::CRYSTALS_IN_TOWER_ETA);
0316           int ref_iPhi = (iLinkPhi * p2eg::CRYSTALS_IN_TOWER_PHI);
0317
0318           // In the link, get the crystals (5x5 in each link)
0319           for (int iEta = 0; iEta < p2eg::CRYSTALS_IN_TOWER_ETA; iEta++) {
0320             for (int iPhi = 0; iPhi < p2eg::CRYSTALS_IN_TOWER_PHI; iPhi++) {
0321               // Et as unsigned int
0322               ap_uint<10> uEnergy = myLink.getCrystalE(iEta, iPhi);
0323
0324               // Fill the 'temporary' array with a crystal object
0325               temporary[ref_iEta + iEta][ref_iPhi + iPhi] = p2eg::crystal(uEnergy);
0326             }
0327           }  // end of loop over crystals
0328
0329           // HCAL tower ET
0330           p2eg::towerHCAL& myTower = myTowers.getTowerHCAL(iLinkEta, iLinkPhi);
0331           towerEtHCAL[(iLinkEta * p2eg::TOWER_IN_PHI) + iLinkPhi] = myTower.getEt();
0332         }
0333       }
0334
0335       // Iteratively find four clusters and remove them from 'temporary' as we go, and fill cluster_list
0336       for (int c = 0; c < p2eg::N_CLUSTERS_PER_REGION; c++) {
0337         p2eg::Cluster newCluster = p2eg::getClusterFromRegion3x4(temporary);  // remove cluster from 'temporary'
0338         newCluster.setRegionIdx(idxRegion);                                   // add the region number
0339         if (newCluster.clusterEnergy() > 0) {
0340           // do not push back 0-energy clusters
0341           cluster_list[cc].push_back(newCluster);
0342         }
0343       }
0344
0345       // Create towers using remaining ECAL energy, and the HCAL towers were already calculated in towersEtHCAL[12]
0346       ap_uint<12> towerEtECAL[12];
0347       p2eg::getECALTowersEt(temporary, towerEtECAL);
0348
0349       // Fill towerHCALCard and towerECALCard arrays
0350       for (int i = 0; i < 12; i++) {
0351         // Get the tower's indices in a (17x4) card
0352         int iEta = (idxRegion * p2eg::TOWER_IN_ETA) + (i / p2eg::TOWER_IN_PHI);
0353         int iPhi = (i % p2eg::TOWER_IN_PHI);
0354
0355         // If the region number is 5 (i.e. only 2x4 in towers, skip the third row) N_REGIONS_PER_CARD = 6. i.e. we do not want to consider
0356         // i = 8, 9, 10, 11
0357         if ((idxRegion == (p2eg::N_REGIONS_PER_CARD - 1)) && (i > 7)) {
0358           continue;
0359         }
0360         towerHCALCard[iEta][iPhi][cc] =
0361             p2eg::tower_t(towerEtHCAL[i], 0);  // p2eg::tower_t initializer takes an ap-uint<12> for the energy
0362         towerECALCard[iEta][iPhi][cc] = p2eg::tower_t(towerEtECAL[i], 0);
0363       }
0364     }
0365
0366     //-------------------------------------------//
0367     // Stitching across ECAL regions             //
0368     //-------------------------------------------//
0369     const int nRegionBoundariesEta = (p2eg::N_REGIONS_PER_CARD - 1);  // 6 regions -> 5 boundaries to check
0370     // Upper and lower boundaries respectively, to check for stitching along
0371     int towerEtaBoundaries[nRegionBoundariesEta][2] = {{15, 14}, {12, 11}, {9, 8}, {6, 5}, {3, 2}};
0372
0373     for (int iBound = 0; iBound < nRegionBoundariesEta; iBound++) {
0374       p2eg::stitchClusterOverRegionBoundary(
0375           cluster_list[cc], towerEtaBoundaries[iBound][0], towerEtaBoundaries[iBound][1], cc);
0376     }
0377
0378     //--------------------------------------------------------------------------------//
0379     // Sort the clusters, take the 8 with highest pT, and return extras to tower
0380     //--------------------------------------------------------------------------------//
0381     if (!cluster_list[cc].empty()) {
0382       std::sort(cluster_list[cc].begin(), cluster_list[cc].end(), p2eg::compareClusterET);
0383
0384       // If there are more than eight clusters, return the unused energy to the towers
0385       for (unsigned int kk = p2eg::n_clusters_4link; kk < cluster_list[cc].size(); ++kk) {
0386         p2eg::Cluster cExtra = cluster_list[cc][kk];
0387         if (cExtra.clusterEnergy() > 0) {
0388           // Increment tower ET
0389           // Get tower (eta, phi) (up to (17, 4)) in the RCT card
0390           int whichTowerEtaInCard = ((cExtra.region() * p2eg::TOWER_IN_ETA) + cExtra.towerEta());
0391           int whichTowerPhiInCard = cExtra.towerPhi();
0392           ap_uint<12> oldTowerEt = towerECALCard[whichTowerEtaInCard][whichTowerPhiInCard][cc].et();
0393           ap_uint<12> newTowerEt = (oldTowerEt + cExtra.clusterEnergy());
0394           ap_uint<4> hoe = towerECALCard[whichTowerEtaInCard][whichTowerPhiInCard][cc].hoe();
0395           towerECALCard[whichTowerEtaInCard][whichTowerPhiInCard][cc] = p2eg::tower_t(newTowerEt, hoe);
0396         }
0397       }
0398
0399       // Save up to eight clusters: loop over cluster_list
0400       for (unsigned int kk = 0; kk < cluster_list[cc].size(); ++kk) {
0401         if (kk >= p2eg::n_clusters_4link)
0402           continue;
0403         if (cluster_list[cc][kk].clusterEnergy() > 0) {
0404           cluster_list_merged[cc].push_back(cluster_list[cc][kk]);
0405         }
0406       }
0407     }
0408
0409     //-------------------------------------------//
0410     // Calibrate clusters
0411     //-------------------------------------------//
0412     for (auto& c : cluster_list_merged[cc]) {
0413       float realEta = c.realEta(cc);
0414       c.calib = calib_(c.getPt(), std::abs(realEta));
0415       c.applyCalibration(c.calib);
0416     }
0417
0418     //-------------------------------------------//
0419     // Cluster shower shape flags
0420     //-------------------------------------------//
0421     for (auto& c : cluster_list_merged[cc]) {
0422       c.is_ss = p2eg::passes_ss(c.getPt(), c.realEta(cc), (c.getEt2x5() / c.getEt5x5()));
0423       c.is_looseTkss = p2eg::passes_looseTkss(c.getPt(), c.realEta(cc), (c.getEt2x5() / c.getEt5x5()));
0424     }
0425
0426     //-------------------------------------------//
0427     // Calibrate towers
0428     //-------------------------------------------//
0429     for (int ii = 0; ii < p2eg::n_towers_cardEta; ++ii) {    // 17 towers per card in eta
0430       for (int jj = 0; jj < p2eg::n_towers_cardPhi; ++jj) {  // 4 towers per card in phi
0431         float tRealEta = p2eg::getTowerEta_fromAbsID(
0432             p2eg::getAbsID_iEta_fromFirmwareCardTowerLink(cc, ii, jj));  // real eta of center of tower
0433         double tCalib = calib_(0, tRealEta);                             // calibration factor
0434         towerECALCard[ii][jj][cc].applyCalibration(tCalib);
0435       }
0436     }
0437
0438     //-------------------------------------------//
0439     // Calculate tower HoE
0440     //-------------------------------------------//
0441     for (int ii = 0; ii < p2eg::n_towers_cardEta; ++ii) {    // 17 towers per card in eta
0442       for (int jj = 0; jj < p2eg::n_towers_cardPhi; ++jj) {  // 4 towers per card in phi
0443         ap_uint<12> ecalEt = towerECALCard[ii][jj][cc].et();
0444         ap_uint<12> hcalEt = towerHCALCard[ii][jj][cc].et();
0445         towerECALCard[ii][jj][cc].addHoverEToTower(ecalEt, hcalEt);
0446       }
0447     }
0448
0449     //-----------------------------------------------------------//
0450     // Produce output RCT collections for event display and analyzer
0451     //-----------------------------------------------------------//
0452     for (auto& c : cluster_list_merged[cc]) {
0453       reco::Candidate::PolarLorentzVector p4calibrated(c.getPt(), c.realEta(cc), c.realPhi(cc), 0.);
0454
0455       // Constructor definition at: https://github.com/cms-l1t-offline/cmssw/blob/l1t-phase2-v3.3.11/DataFormats/L1TCalorimeterPhase2/interface/CaloCrystalCluster.h#L34
0456       l1tp2::CaloCrystalCluster cluster(p4calibrated,
0457                                         c.getPt(),           // use float
0458                                         0,                   // float h over e
0459                                         0,                   // float iso
0460                                         0,                   // DetId seedCrystal
0461                                         0,                   // puCorrPt
0462                                         c.getBrems(),        // 0, 1, or 2 (as computed in firmware)
0463                                         0,                   // et2x2 (not calculated)
0464                                         c.getEt2x5(),        // et2x5 (as computed in firmware, save float)
0465                                         0,                   // et3x5 (not calculated)
0466                                         c.getEt5x5(),        // et5x5 (as computed in firmware, save float)
0467                                         c.getIsSS(),         // standalone WP
0468                                         c.getIsSS(),         // electronWP98
0469                                         false,               // photonWP80
0470                                         c.getIsSS(),         // electronWP90
0471                                         c.getIsLooseTkss(),  // looseL1TkMatchWP
0472                                         c.getIsSS()          // stage2effMatch
0473       );
0474
0475       std::map<std::string, float> params;
0476       params["standaloneWP_showerShape"] = c.getIsSS();
0477       params["trkMatchWP_showerShape"] = c.getIsLooseTkss();
0478       cluster.setExperimentalParams(params);
0479
0480       L1EGXtalClusters->push_back(cluster);
0481     }
0482     // Output tower collections
0483     for (int ii = 0; ii < p2eg::n_towers_cardEta; ++ii) {    // 17 towers per card in eta
0484       for (int jj = 0; jj < p2eg::n_towers_cardPhi; ++jj) {  // 4 towers per card in phi
0485
0486         l1tp2::CaloTower l1CaloTower;
0487         l1CaloTower.setEcalTowerEt(towerECALCard[ii][jj][cc].et() * p2eg::ECAL_LSB);
0488         l1CaloTower.setHcalTowerEt(towerHCALCard[ii][jj][cc].et() * p2eg::HCAL_LSB);
0489         int absToweriEta = p2eg::getAbsID_iEta_fromFirmwareCardTowerLink(cc, ii, jj);
0490         int absToweriPhi = p2eg::getAbsID_iPhi_fromFirmwareCardTowerLink(cc, ii, jj);
0491         l1CaloTower.setTowerIEta(absToweriEta);
0492         l1CaloTower.setTowerIPhi(absToweriPhi);
0493         l1CaloTower.setTowerEta(p2eg::getTowerEta_fromAbsID(absToweriEta));
0494         l1CaloTower.setTowerPhi(p2eg::getTowerPhi_fromAbsID(absToweriPhi));
0495
0496         L1CaloTowers->push_back(l1CaloTower);
0497       }
0498     }
0499   }  // end of loop over cards
0500
0501   iEvent.put(std::move(L1EGXtalClusters), "RCTClusters");
0502   iEvent.put(std::move(L1CaloTowers), "RCTTowers");
0503
0504   //*******************************************************************
0505   // Do GCT geometry for inputs
0506   //*******************************************************************
0507
0508   // Loop over GCT cards (three of them)
0509   p2eg::GCTcard_t gctCards[p2eg::N_GCTCARDS];
0510   p2eg::GCTtoCorr_t gctToCorr[p2eg::N_GCTCARDS];
0511
0512   // Initialize the cards (requires towerECALCard, towerHCALCard arrays, and cluster_list_merged)
0513   for (unsigned int gcc = 0; gcc < p2eg::N_GCTCARDS; gcc++) {
0514     // Each GCT card encompasses 16 RCT cards, listed in GCTcardtoRCTcardnumber[3][16]. i goes from 0 to <16
0515     for (int i = 0; i < (p2eg::N_RCTCARDS_PHI * 2); i++) {
0516       unsigned int rcc = p2eg::GCTcardtoRCTcardnumber[gcc][i];
0517
0518       // Positive eta? Fist row is in positive eta
0519       bool isPositiveEta = (i < p2eg::N_RCTCARDS_PHI);
0520
0521       // Sum tower ECAL and HCAL energies: 17 towers per link
0522       for (int iTower = 0; iTower < p2eg::N_GCTTOWERS_FIBER; iTower++) {
0523         // 4 links per card
0524         for (int iLink = 0; iLink < p2eg::n_links_card; iLink++) {
0525           p2eg::tower_t t0_ecal = towerECALCard[iTower][iLink][rcc];
0526           p2eg::tower_t t0_hcal = towerHCALCard[iTower][iLink][rcc];
0527           p2eg::RCTtower_t t;
0528           t.et = t0_ecal.et() + t0_hcal.et();
0529           t.hoe = t0_ecal.hoe();
0530           // Not needed for GCT firmware but will be written into GCT CMSSW outputs : 12 bits each
0531           t.ecalEt = t0_ecal.et();
0532           t.hcalEt = t0_hcal.et();
0533
0534           if (isPositiveEta) {
0535             gctCards[gcc].RCTcardEtaPos[i % p2eg::N_RCTCARDS_PHI].RCTtoGCTfiber[iLink].RCTtowers[iTower] = t;
0536           } else {
0537             gctCards[gcc].RCTcardEtaNeg[i % p2eg::N_RCTCARDS_PHI].RCTtoGCTfiber[iLink].RCTtowers[iTower] = t;
0538           }
0539         }
0540       }
0541
0542       // Distribute at most 8 RCT clusters across four links: convert to GCT coordinates
0543       for (size_t iCluster = 0; (iCluster < cluster_list_merged[rcc].size()) &&
0544                                 (iCluster < (p2eg::N_RCTGCT_FIBERS * p2eg::N_RCTCLUSTERS_FIBER));
0545            iCluster++) {
0546         p2eg::Cluster c0 = cluster_list_merged[rcc][iCluster];
0547         p2eg::RCTcluster_t c;
0548         c.et = c0.clusterEnergy();
0549
0550         // tower Eta: c0.towerEta() refers to the tower iEta INSIDE the region, so we need to convert to tower iEta inside the card
0551         c.towEta = (c0.region() * p2eg::TOWER_IN_ETA) + c0.towerEta();
0552         c.towPhi = c0.towerPhi();
0553         c.crEta = c0.clusterEta();
0554         c.crPhi = c0.clusterPhi();
0555         c.et5x5 = c0.uint_et5x5();
0556         c.et2x5 = c0.uint_et2x5();
0557         c.is_ss = c0.getIsSS();
0558         c.is_looseTkss = c0.getIsLooseTkss();
0559         c.is_iso = c0.getIsIso();
0560         c.is_looseTkiso = c0.getIsLooseTkIso();
0561         c.brems = c0.getBrems();
0562         c.nGCTCard = gcc;  // store gct card index as well
0563         unsigned int iIdxInGCT = i % p2eg::N_RCTCARDS_PHI;
0564         unsigned int iLinkC = iCluster % p2eg::N_RCTGCT_FIBERS;
0565         unsigned int iPosC = iCluster / p2eg::N_RCTGCT_FIBERS;
0566
0567         if (isPositiveEta) {
0568           gctCards[gcc].RCTcardEtaPos[iIdxInGCT].RCTtoGCTfiber[iLinkC].RCTclusters[iPosC] = c;
0569         } else {
0570           gctCards[gcc].RCTcardEtaNeg[iIdxInGCT].RCTtoGCTfiber[iLinkC].RCTclusters[iPosC] = c;
0571         }
0572       }
0573       // If there were fewer than eight clusters, make sure the remaining fiber clusters are zero'd out.
0574       for (size_t iZeroCluster = cluster_list_merged[rcc].size();
0575            iZeroCluster < (p2eg::N_RCTGCT_FIBERS * p2eg::N_RCTCLUSTERS_FIBER);
0576            iZeroCluster++) {
0577         unsigned int iIdxInGCT = i % p2eg::N_RCTCARDS_PHI;
0578         unsigned int iLinkC = iZeroCluster % p2eg::N_RCTGCT_FIBERS;
0579         unsigned int iPosC = iZeroCluster / p2eg::N_RCTGCT_FIBERS;
0580
0581         p2eg::RCTcluster_t cZero;
0582         cZero.et = 0;
0583         cZero.towEta = 0;
0584         cZero.towPhi = 0;
0585         cZero.crEta = 0;
0586         cZero.crPhi = 0;
0587         cZero.et5x5 = 0;
0588         cZero.et2x5 = 0;
0589         cZero.is_ss = false;
0590         cZero.is_looseTkss = false;
0591         cZero.is_iso = false;
0592         cZero.is_looseTkiso = false;
0593         cZero.nGCTCard = gcc;  // store gct card index as well
0594         if (isPositiveEta) {
0595           gctCards[gcc].RCTcardEtaPos[iIdxInGCT].RCTtoGCTfiber[iLinkC].RCTclusters[iPosC] = cZero;
0596         } else {
0597           gctCards[gcc].RCTcardEtaNeg[iIdxInGCT].RCTtoGCTfiber[iLinkC].RCTclusters[iPosC] = cZero;
0598         }
0599       }
0600     }
0601   }  // end of loop over initializing GCT cards
0602
0603   //----------------------------------------------------
0604   // Output collections for the GCT emulator
0605   //----------------------------------------------------
0606   auto L1GCTClusters = std::make_unique<l1tp2::CaloCrystalClusterCollection>();
0607   auto L1GCTTowers = std::make_unique<l1tp2::CaloTowerCollection>();
0608   auto L1GCTFullTowers = std::make_unique<l1tp2::CaloTowerCollection>();
0609   auto L1GCTEGammas = std::make_unique<l1t::EGammaBxCollection>();
0610   auto L1DigitizedClusterCorrelator = std::make_unique<l1tp2::DigitizedClusterCorrelatorCollection>();
0611   auto L1DigitizedTowerCorrelator = std::make_unique<l1tp2::DigitizedTowerCorrelatorCollection>();
0612   auto L1DigitizedClusterGT = std::make_unique<l1tp2::DigitizedClusterGTCollection>();
0613
0614   //----------------------------------------------------
0615   // Apply the GCT firmware code to each GCT card
0616   //----------------------------------------------------
0617
0618   for (unsigned int gcc = 0; gcc < p2eg::N_GCTCARDS; gcc++) {
0619     p2eg::algo_top(gctCards[gcc],
0620                    gctToCorr[gcc],
0621                    gcc,
0622                    L1GCTClusters,
0623                    L1GCTTowers,
0624                    L1GCTFullTowers,
0625                    L1GCTEGammas,
0626                    L1DigitizedClusterCorrelator,
0627                    L1DigitizedTowerCorrelator,
0628                    L1DigitizedClusterGT,
0629                    calib_);
0630   }
0631
0632   iEvent.put(std::move(L1GCTClusters), "GCTClusters");
0633   iEvent.put(std::move(L1GCTTowers), "GCTTowers");
0634   iEvent.put(std::move(L1GCTFullTowers), "GCTFullTowers");
0635   iEvent.put(std::move(L1GCTEGammas), "GCTEGammas");
0636   iEvent.put(std::move(L1DigitizedClusterCorrelator), "GCTDigitizedClusterToCorrelator");
0637   iEvent.put(std::move(L1DigitizedTowerCorrelator), "GCTDigitizedTowerToCorrelator");
0638   iEvent.put(std::move(L1DigitizedClusterGT), "GCTDigitizedClusterToGT");
0639 }
0640
0641 //////////////////////////////////////////////////////////////////////////
0642
0643 void Phase2L1CaloEGammaEmulator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0644   // l1tPhase2L1CaloEGammaEmulator
0645   edm::ParameterSetDescription desc;
0646   desc.add<edm::InputTag>("ecalTPEB", edm::InputTag("simEcalEBTriggerPrimitiveDigis"));
0647   desc.add<edm::InputTag>("hcalTP", edm::InputTag("simHcalTriggerPrimitiveDigis"));
0648   {
0649     edm::ParameterSetDescription psd0;
0650     psd0.add<std::vector<double>>("etaBins",
0651                                   {
0652                                       0.087,
0653                                       0.174,
0654                                       0.261,
0655                                       0.348,
0656                                       0.435,
0657                                       0.522,
0658                                       0.609,
0659                                       0.696,
0660                                       0.783,
0661                                       0.87,
0662                                       0.957,
0663                                       1.044,
0664                                       1.131,
0665                                       1.218,
0666                                       1.305,
0667                                       1.392,
0668                                       1.479,
0669                                   });
0670     psd0.add<std::vector<double>>("ptBins",
0671                                   {
0672                                       12,
0673                                       20,
0674                                       30,
0675                                       40,
0676                                       55,
0677                                       90,
0678                                       1000000.0,
0679                                   });
0680     psd0.add<std::vector<double>>("scale",
0681                                   {
0682                                       1.298,  1.287,
0683                                       1.309,  1.298,
0684                                       1.309,  1.309,
0685                                       1.309,  1.298,
0686                                       1.309,  1.298,
0687                                       1.309,  1.309,
0688                                       1.309,  1.32,
0689                                       1.309,  1.32,
0690                                       1.309,  1.1742,
0691                                       1.1639, 1.1639,
0692                                       1.1639, 1.1639,
0693                                       1.1639, 1.1639,
0694                                       1.1742, 1.1742,
0695                                       1.1639, 1.1639,
0696                                       1.1742, 1.1639,
0697                                       1.1639, 1.1742,
0698                                       1.1742, 1.1536000000000002,
0699                                       1.11,   1.11,
0700                                       1.11,   1.11,
0701                                       1.11,   1.11,
0702                                       1.11,   1.11,
0703                                       1.11,   1.11,
0704                                       1.11,   1.11,
0705                                       1.11,   1.11,
0706                                       1.11,   1.11,
0707                                       1.1,    1.09,
0708                                       1.09,   1.09,
0709                                       1.09,   1.09,
0710                                       1.09,   1.09,
0711                                       1.09,   1.09,
0712                                       1.09,   1.09,
0713                                       1.09,   1.09,
0714                                       1.09,   1.09,
0715                                       1.09,   1.09,
0716                                       1.07,   1.07,
0717                                       1.07,   1.07,
0718                                       1.07,   1.07,
0719                                       1.07,   1.08,
0720                                       1.07,   1.07,
0721                                       1.08,   1.08,
0722                                       1.07,   1.08,
0723                                       1.08,   1.08,
0724                                       1.08,   1.06,
0725                                       1.06,   1.06,
0726                                       1.06,   1.05,
0727                                       1.05,   1.06,
0728                                       1.06,   1.06,
0729                                       1.06,   1.06,
0730                                       1.06,   1.06,
0731                                       1.06,   1.06,
0732                                       1.06,   1.06,
0733                                       1.04,   1.04,
0734                                       1.04,   1.04,
0735                                       1.05,   1.04,
0736                                       1.05,   1.05,
0737                                       1.05,   1.05,
0738                                       1.05,   1.05,
0739                                       1.05,   1.05,
0740                                       1.05,   1.05,
0741                                       1.05,
0742                                   });
0743     desc.add<edm::ParameterSetDescription>("calib", psd0);
0744   }
0745   descriptions.addWithDefaultLabel(desc);
0746 }
0747
0748 //define this as a plug-in
0749 DEFINE_FWK_MODULE(Phase2L1CaloEGammaEmulator);