Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:54:43 

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>("RCT");
0096   produces<l1tp2::CaloCrystalClusterCollection>("GCT");
0097   produces<l1tp2::CaloTowerCollection>("RCT");
0098   produces<l1tp2::CaloTowerCollection>("GCT");
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 divide by 8
0138       float et = hit.encodedEt() * p2eg::ECAL_LSB;
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((ap_uint<10>)hit.encodedEt());  // also save the uint Et
0151       ehit.setPt();
0152       ecalhits.push_back(ehit);
0153     }
0154   }
0155 
0156   //***************************************************//
0157   // Get the HCAL hits
0158   //***************************************************//
0159   std::vector<p2eg::SimpleCaloHit> hcalhits;
0160   edm::Handle<edm::SortedCollection<HcalTriggerPrimitiveDigi>> hbhecoll;
0161   iEvent.getByToken(hcalTPToken_, hbhecoll);
0162 
0163   for (const auto& hit : *hbhecoll.product()) {
0164     float et = decoder.hcaletValue(hit.id(), hit.t0());
0165     ap_uint<10> encodedEt = hit.t0().compressedEt();
0166     // same thing as SOI_compressedEt() in HcalTriggerPrimitiveDigi.h///
0167     if (et <= 0)
0168       continue;
0169 
0170     if (!(hcTopology_->validHT(hit.id()))) {
0171       LogError("Phase2L1CaloEGammaEmulator")
0172           << " -- Hcal hit DetID not present in HCAL Geom: " << hit.id() << std::endl;
0173       throw cms::Exception("Phase2L1CaloEGammaEmulator");
0174       continue;
0175     }
0176     const std::vector<HcalDetId>& hcId = theTrigTowerGeometry.detIds(hit.id());
0177     if (hcId.empty()) {
0178       LogError("Phase2L1CaloEGammaEmulator") << "Cannot find any HCalDetId corresponding to " << hit.id() << std::endl;
0179       throw cms::Exception("Phase2L1CaloEGammaEmulator");
0180       continue;
0181     }
0182     if (hcId[0].subdetId() > 1) {
0183       continue;
0184     }
0185     GlobalVector hcal_tp_position = GlobalVector(0., 0., 0.);
0186     for (const auto& hcId_i : hcId) {
0187       if (hcId_i.subdetId() > 1) {
0188         continue;
0189       }
0190       // get the first HCAL TP/ cell
0191       auto cell = hbGeometry->getGeometry(hcId_i);
0192       if (cell == nullptr) {
0193         continue;
0194       }
0195       GlobalVector tmpVector = GlobalVector(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z());
0196       hcal_tp_position = tmpVector;
0197 
0198       break;
0199     }
0200     p2eg::SimpleCaloHit hhit;
0201     hhit.setId(hit.id());
0202     hhit.setIdHcal(hit.id());
0203     hhit.setPosition(hcal_tp_position);
0204     hhit.setEnergy(et);
0205     hhit.setPt();
0206     hhit.setEt_uint(encodedEt);
0207     hcalhits.push_back(hhit);
0208   }
0209 
0210   //***************************************************//
0211   // Initialize necessary arrays for tower and clusters
0212   //***************************************************//
0213 
0214   // L1 Outputs definition: Arrays that use firmware convention for indexing
0215   p2eg::tower_t towerHCALCard
0216       [p2eg::n_towers_cardEta][p2eg::n_towers_cardPhi]
0217       [p2eg::n_towers_halfPhi];  // 17x4x36 array (not to be confused with the 12x1 array of ap_uints, towerEtHCAL
0218   p2eg::tower_t towerECALCard[p2eg::n_towers_cardEta][p2eg::n_towers_cardPhi][p2eg::n_towers_halfPhi];
0219   // There is one vector of clusters per card (up to 12 clusters before stitching across ECAL regions)
0220   std::vector<p2eg::Cluster> cluster_list[p2eg::n_towers_halfPhi];
0221   // After merging/stitching the clusters, we only take the 8 highest pt per card
0222   std::vector<p2eg::Cluster> cluster_list_merged[p2eg::n_towers_halfPhi];
0223 
0224   //***************************************************//
0225   // Fill RCT ECAL regions with ECAL hits
0226   //***************************************************//
0227   for (int cc = 0; cc < p2eg::n_towers_halfPhi; ++cc) {  // Loop over 36 L1 cards
0228 
0229     p2eg::card rctCard;
0230     rctCard.setIdx(cc);
0231 
0232     for (const auto& hit : ecalhits) {
0233       // Check if the hit is in cards 0-35
0234       if (hit.isInCard(cc)) {
0235         // Get the crystal eta and phi, relative to the bottom left corner of the card
0236         // (0 up to 17*5, 0 up to 4*5)
0237         int local_iEta = hit.crystalLocaliEta(cc);
0238         int local_iPhi = hit.crystalLocaliPhi(cc);
0239 
0240         // Region number (0-5) depends only on the crystal iEta in the card
0241         int regionNumber = p2eg::getRegionNumber(local_iEta);
0242 
0243         // Tower eta and phi index inside the card (17x4)
0244         int inCard_tower_iEta = local_iEta / p2eg::CRYSTALS_IN_TOWER_ETA;
0245         int inCard_tower_iPhi = local_iPhi / p2eg::CRYSTALS_IN_TOWER_PHI;
0246 
0247         // Tower eta and phi index inside the region (3x4)
0248         int inRegion_tower_iEta = inCard_tower_iEta % p2eg::TOWER_IN_ETA;
0249         int inRegion_tower_iPhi = inCard_tower_iPhi % p2eg::TOWER_IN_PHI;
0250 
0251         // Crystal eta and phi index inside the 3x4 region (15x20)
0252         int inRegion_crystal_iEta = local_iEta % (p2eg::TOWER_IN_ETA * p2eg::CRYSTALS_IN_TOWER_ETA);
0253         int inRegion_crystal_iPhi = local_iPhi;
0254 
0255         // Crystal eta and phi index inside the tower (5x5)
0256         int inLink_crystal_iEta = (inRegion_crystal_iEta % p2eg::CRYSTALS_IN_TOWER_ETA);
0257         int inLink_crystal_iPhi = (inRegion_crystal_iPhi % p2eg::CRYSTALS_IN_TOWER_PHI);
0258 
0259         // Add the crystal energy to the rctCard
0260         p2eg::region3x4& myRegion = rctCard.getRegion3x4(regionNumber);
0261         p2eg::linkECAL& myLink = myRegion.getLinkECAL(inRegion_tower_iEta, inRegion_tower_iPhi);
0262         myLink.addCrystalE(inLink_crystal_iEta, inLink_crystal_iPhi, hit.et_uint());
0263       }
0264     }
0265 
0266     //***************************************************//
0267     // Build RCT towers from HCAL hits
0268     //***************************************************//
0269     for (const auto& hit : hcalhits) {
0270       if (hit.isInCard(cc) && hit.pt() > 0) {
0271         // Get crystal eta and phi, relative to the bottom left corner of the card
0272         // (0 up to 17*5, 0 up to 4*5)
0273         int local_iEta = hit.crystalLocaliEta(cc);
0274         int local_iPhi = hit.crystalLocaliPhi(cc);
0275 
0276         // Region (0-5) the hit falls into
0277         int regionNumber = p2eg::getRegionNumber(local_iEta);
0278 
0279         // Tower eta and phi index inside the card (17x4)
0280         int inCard_tower_iEta = int(local_iEta / p2eg::CRYSTALS_IN_TOWER_ETA);
0281         int inCard_tower_iPhi = int(local_iPhi / p2eg::CRYSTALS_IN_TOWER_PHI);
0282 
0283         // Tower eta and phi index inside the region (3x4)
0284         int inRegion_tower_iEta = inCard_tower_iEta % p2eg::TOWER_IN_ETA;
0285         int inRegion_tower_iPhi = inCard_tower_iPhi % p2eg::TOWER_IN_PHI;
0286 
0287         // Access the right HCAL region and tower and increment the ET
0288         p2eg::towers3x4& myTowers3x4 = rctCard.getTowers3x4(regionNumber);
0289         p2eg::towerHCAL& myTower = myTowers3x4.getTowerHCAL(inRegion_tower_iEta, inRegion_tower_iPhi);
0290         myTower.addEt(hit.et_uint());
0291       }
0292     }
0293 
0294     //***************************************************//
0295     // Make clusters in each ECAL region independently
0296     //***************************************************//
0297     for (int idxRegion = 0; idxRegion < p2eg::N_REGIONS_PER_CARD; idxRegion++) {
0298       // ECAL crystals array
0299       p2eg::crystal temporary[p2eg::CRYSTAL_IN_ETA][p2eg::CRYSTAL_IN_PHI];
0300       // HCAL towers in 3x4 region
0301       ap_uint<12> towerEtHCAL[p2eg::TOWER_IN_ETA * p2eg::TOWER_IN_PHI];
0302 
0303       p2eg::region3x4& myRegion = rctCard.getRegion3x4(idxRegion);
0304       p2eg::towers3x4& myTowers = rctCard.getTowers3x4(idxRegion);
0305 
0306       // In each 3x4 region, loop through the links (one link = one tower)
0307       for (int iLinkEta = 0; iLinkEta < p2eg::TOWER_IN_ETA; iLinkEta++) {
0308         for (int iLinkPhi = 0; iLinkPhi < p2eg::TOWER_IN_PHI; iLinkPhi++) {
0309           // Get the ECAL link (one link per tower)
0310           p2eg::linkECAL& myLink = myRegion.getLinkECAL(iLinkEta, iLinkPhi);
0311 
0312           // 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
0313           int ref_iEta = (iLinkEta * p2eg::CRYSTALS_IN_TOWER_ETA);
0314           int ref_iPhi = (iLinkPhi * p2eg::CRYSTALS_IN_TOWER_PHI);
0315 
0316           // In the link, get the crystals (5x5 in each link)
0317           for (int iEta = 0; iEta < p2eg::CRYSTALS_IN_TOWER_ETA; iEta++) {
0318             for (int iPhi = 0; iPhi < p2eg::CRYSTALS_IN_TOWER_PHI; iPhi++) {
0319               // Et as unsigned int
0320               ap_uint<10> uEnergy = myLink.getCrystalE(iEta, iPhi);
0321 
0322               // Fill the 'temporary' array with a crystal object
0323               temporary[ref_iEta + iEta][ref_iPhi + iPhi] = p2eg::crystal(uEnergy);
0324             }
0325           }  // end of loop over crystals
0326 
0327           // HCAL tower ET
0328           p2eg::towerHCAL& myTower = myTowers.getTowerHCAL(iLinkEta, iLinkPhi);
0329           towerEtHCAL[(iLinkEta * p2eg::TOWER_IN_PHI) + iLinkPhi] = myTower.getEt();
0330         }
0331       }
0332 
0333       // Iteratively find four clusters and remove them from 'temporary' as we go, and fill cluster_list
0334       for (int c = 0; c < p2eg::N_CLUSTERS_PER_REGION; c++) {
0335         p2eg::Cluster newCluster = p2eg::getClusterFromRegion3x4(temporary);  // remove cluster from 'temporary'
0336         newCluster.setRegionIdx(idxRegion);                                   // add the region number
0337         if (newCluster.clusterEnergy() > 0) {
0338           // do not push back 0-energy clusters
0339           cluster_list[cc].push_back(newCluster);
0340         }
0341       }
0342 
0343       // Create towers using remaining ECAL energy, and the HCAL towers were already calculated in towersEtHCAL[12]
0344       ap_uint<12> towerEtECAL[12];
0345       p2eg::getECALTowersEt(temporary, towerEtECAL);
0346 
0347       // Fill towerHCALCard and towerECALCard arrays
0348       for (int i = 0; i < 12; i++) {
0349         // Get the tower's indices in a (17x4) card
0350         int iEta = (idxRegion * p2eg::TOWER_IN_ETA) + (i / p2eg::TOWER_IN_PHI);
0351         int iPhi = (i % p2eg::TOWER_IN_PHI);
0352 
0353         // 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
0354         // i = 8, 9, 10, 11
0355         if ((idxRegion == (p2eg::N_REGIONS_PER_CARD - 1)) && (i > 7)) {
0356           continue;
0357         }
0358         towerHCALCard[iEta][iPhi][cc] =
0359             p2eg::tower_t(towerEtHCAL[i], 0);  // p2eg::tower_t initializer takes an ap-uint<12> for the energy
0360         towerECALCard[iEta][iPhi][cc] = p2eg::tower_t(towerEtECAL[i], 0);
0361       }
0362     }
0363 
0364     //-------------------------------------------//
0365     // Stitching across ECAL regions             //
0366     //-------------------------------------------//
0367     const int nRegionBoundariesEta = (p2eg::N_REGIONS_PER_CARD - 1);  // 6 regions -> 5 boundaries to check
0368     // Upper and lower boundaries respectively, to check for stitching along
0369     int towerEtaBoundaries[nRegionBoundariesEta][2] = {{15, 14}, {12, 11}, {9, 8}, {6, 5}, {3, 2}};
0370 
0371     for (int iBound = 0; iBound < nRegionBoundariesEta; iBound++) {
0372       p2eg::stitchClusterOverRegionBoundary(
0373           cluster_list[cc], towerEtaBoundaries[iBound][0], towerEtaBoundaries[iBound][1], cc);
0374     }
0375 
0376     //--------------------------------------------------------------------------------//
0377     // Sort the clusters, take the 8 with highest pT, and return extras to tower
0378     //--------------------------------------------------------------------------------//
0379     if (!cluster_list[cc].empty()) {
0380       std::sort(cluster_list[cc].begin(), cluster_list[cc].end(), p2eg::compareClusterET);
0381 
0382       // If there are more than eight clusters, return the unused energy to the towers
0383       for (unsigned int kk = p2eg::n_clusters_4link; kk < cluster_list[cc].size(); ++kk) {
0384         p2eg::Cluster cExtra = cluster_list[cc][kk];
0385         if (cExtra.clusterEnergy() > 0) {
0386           // Increment tower ET
0387           // Get tower (eta, phi) (up to (17, 4)) in the RCT card
0388           int whichTowerEtaInCard = ((cExtra.region() * p2eg::TOWER_IN_ETA) + cExtra.towerEta());
0389           int whichTowerPhiInCard = cExtra.towerPhi();
0390           ap_uint<12> oldTowerEt = towerECALCard[whichTowerEtaInCard][whichTowerPhiInCard][cc].et();
0391           ap_uint<12> newTowerEt = (oldTowerEt + cExtra.clusterEnergy());
0392           ap_uint<4> hoe = towerECALCard[whichTowerEtaInCard][whichTowerPhiInCard][cc].hoe();
0393           towerECALCard[whichTowerEtaInCard][whichTowerPhiInCard][cc] = p2eg::tower_t(newTowerEt, hoe);
0394         }
0395       }
0396 
0397       // Save up to eight clusters: loop over cluster_list
0398       for (unsigned int kk = 0; kk < cluster_list[cc].size(); ++kk) {
0399         if (kk >= p2eg::n_clusters_4link)
0400           continue;
0401         if (cluster_list[cc][kk].clusterEnergy() > 0) {
0402           cluster_list_merged[cc].push_back(cluster_list[cc][kk]);
0403         }
0404       }
0405     }
0406 
0407     //-------------------------------------------//
0408     // Calibrate clusters
0409     //-------------------------------------------//
0410     for (auto& c : cluster_list_merged[cc]) {
0411       float realEta = c.realEta(cc);
0412       c.calib = calib_(c.getPt(), std::abs(realEta));
0413       c.applyCalibration(c.calib);
0414     }
0415 
0416     //-------------------------------------------//
0417     // Cluster shower shape flags
0418     //-------------------------------------------//
0419     for (auto& c : cluster_list_merged[cc]) {
0420       c.is_ss = p2eg::passes_ss(c.getPt(), (c.getEt2x5() / c.getEt5x5()));
0421       c.is_looseTkss = p2eg::passes_looseTkss(c.getPt(), (c.getEt2x5() / c.getEt5x5()));
0422     }
0423 
0424     //-------------------------------------------//
0425     // Calibrate towers
0426     //-------------------------------------------//
0427     for (int ii = 0; ii < p2eg::n_towers_cardEta; ++ii) {    // 17 towers per card in eta
0428       for (int jj = 0; jj < p2eg::n_towers_cardPhi; ++jj) {  // 4 towers per card in phi
0429         float tRealEta = p2eg::getTowerEta_fromAbsID(
0430             p2eg::getAbsID_iEta_fromFirmwareCardTowerLink(cc, ii, jj));  // real eta of center of tower
0431         double tCalib = calib_(0, tRealEta);                             // calibration factor
0432         towerECALCard[ii][jj][cc].applyCalibration(tCalib);
0433       }
0434     }
0435 
0436     //-------------------------------------------//
0437     // Calculate tower HoE
0438     //-------------------------------------------//
0439     for (int ii = 0; ii < p2eg::n_towers_cardEta; ++ii) {    // 17 towers per card in eta
0440       for (int jj = 0; jj < p2eg::n_towers_cardPhi; ++jj) {  // 4 towers per card in phi
0441         ap_uint<12> ecalEt = towerECALCard[ii][jj][cc].et();
0442         ap_uint<12> hcalEt = towerHCALCard[ii][jj][cc].et();
0443         towerECALCard[ii][jj][cc].getHoverE(ecalEt, hcalEt);
0444       }
0445     }
0446 
0447     //-----------------------------------------------------------//
0448     // Produce output RCT collections for event display and analyzer
0449     //-----------------------------------------------------------//
0450     for (auto& c : cluster_list_merged[cc]) {
0451       reco::Candidate::PolarLorentzVector p4calibrated(c.getPt(), c.realEta(cc), c.realPhi(cc), 0.);
0452 
0453       // Constructor definition at: https://github.com/cms-l1t-offline/cmssw/blob/l1t-phase2-v3.3.11/DataFormats/L1TCalorimeterPhase2/interface/CaloCrystalCluster.h#L34
0454       l1tp2::CaloCrystalCluster cluster(p4calibrated,
0455                                         c.getPt(),           // use float
0456                                         0,                   // float h over e
0457                                         0,                   // float iso
0458                                         0,                   // DetId seedCrystal
0459                                         0,                   // puCorrPt
0460                                         c.getBrems(),        // 0, 1, or 2 (as computed in firmware)
0461                                         0,                   // et2x2 (not calculated)
0462                                         c.getEt2x5(),        // et2x5 (as computed in firmware, save float)
0463                                         0,                   // et3x5 (not calculated)
0464                                         c.getEt5x5(),        // et5x5 (as computed in firmware, save float)
0465                                         c.getIsSS(),         // standalone WP
0466                                         c.getIsSS(),         // electronWP98
0467                                         false,               // photonWP80
0468                                         c.getIsSS(),         // electronWP90
0469                                         c.getIsLooseTkss(),  // looseL1TkMatchWP
0470                                         c.getIsSS()          // stage2effMatch
0471       );
0472 
0473       std::map<std::string, float> params;
0474       params["standaloneWP_showerShape"] = c.getIsSS();
0475       params["trkMatchWP_showerShape"] = c.getIsLooseTkss();
0476       cluster.setExperimentalParams(params);
0477 
0478       L1EGXtalClusters->push_back(cluster);
0479     }
0480     // Output tower collections
0481     for (int ii = 0; ii < p2eg::n_towers_cardEta; ++ii) {    // 17 towers per card in eta
0482       for (int jj = 0; jj < p2eg::n_towers_cardPhi; ++jj) {  // 4 towers per card in phi
0483 
0484         l1tp2::CaloTower l1CaloTower;
0485         // Divide by 8.0 to get ET as float (GeV)
0486         l1CaloTower.setEcalTowerEt(towerECALCard[ii][jj][cc].et() * p2eg::ECAL_LSB);
0487         // HCAL TPGs encoded ET: multiply by the LSB (0.5) to convert to GeV
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), "RCT");
0502   iEvent.put(std::move(L1CaloTowers), "RCT");
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() + p2eg::convertHcalETtoEcalET(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), "GCT");
0633   iEvent.put(std::move(L1GCTTowers), "GCT");
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!