Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-09-05 05:04:37

0001 // -*- C++ -*-
0002 //
0003 // Package:    L1Trigger/L1CaloTrigger
0004 // Class:      Phase2L1CaloPFClusterEmulator
0005 //
0006 /**\class Phase2L1CaloPFClusterEmulator Phase2L1CaloPFClusterEmulator.cc L1Trigger/L1CaloTrigger/plugins/Phase2L1CaloPFClusterEmulator.cc
0007 
0008  Description: Creates 3x3 PF clusters from GCTintTowers to be sent to correlator. Follows firmware logic, creates 8 clusters per (2+17+2)x(2+4+2).
0009 
0010  Implementation: To be run together with Phase2L1CaloEGammaEmulator.
0011 
0012 */
0013 //
0014 // Original Author:  Pallabi Das
0015 //         Created:  Wed, 21 Sep 2022 14:54:20 GMT
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/stream/EDProducer.h"
0025 
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/Utilities/interface/StreamID.h"
0031 
0032 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloCrystalCluster.h"
0033 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloTower.h"
0034 #include "DataFormats/L1TCalorimeterPhase2/interface/CaloPFCluster.h"
0035 #include "DataFormats/L1Trigger/interface/EGamma.h"
0036 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0037 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0038 #include "CalibFormats/CaloTPG/interface/CaloTPGTranscoder.h"
0039 #include "CalibFormats/CaloTPG/interface/CaloTPGRecord.h"
0040 #include "L1Trigger/L1TCalorimeter/interface/CaloTools.h"
0041 
0042 #include <ap_int.h>
0043 #include <fstream>
0044 #include <iomanip>
0045 #include <iostream>
0046 #include <cstdio>
0047 #include "L1Trigger/L1CaloTrigger/interface/Phase2L1CaloPFClusterEmulator.h"
0048 
0049 //
0050 // class declaration
0051 //
0052 
0053 class Phase2L1CaloPFClusterEmulator : public edm::stream::EDProducer<> {
0054 public:
0055   explicit Phase2L1CaloPFClusterEmulator(const edm::ParameterSet&);
0056   ~Phase2L1CaloPFClusterEmulator() override = default;
0057 
0058   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0059 
0060 private:
0061   void produce(edm::Event&, const edm::EventSetup&) override;
0062 
0063   // ----------member data ---------------------------
0064   edm::EDGetTokenT<l1tp2::CaloTowerCollection> caloTowerToken_;
0065   edm::EDGetTokenT<HcalTrigPrimDigiCollection> hfToken_;
0066   edm::ESGetToken<CaloTPGTranscoder, CaloTPGRecord> decoderTag_;
0067 };
0068 
0069 //
0070 // constructors and destructor
0071 //
0072 Phase2L1CaloPFClusterEmulator::Phase2L1CaloPFClusterEmulator(const edm::ParameterSet& iConfig)
0073     : caloTowerToken_(consumes<l1tp2::CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("gctFullTowers"))),
0074       hfToken_(consumes<HcalTrigPrimDigiCollection>(iConfig.getParameter<edm::InputTag>("hcalDigis"))),
0075       decoderTag_(esConsumes<CaloTPGTranscoder, CaloTPGRecord>(edm::ESInputTag("", ""))) {
0076   produces<l1tp2::CaloPFClusterCollection>("GCTPFCluster");
0077 }
0078 
0079 //
0080 // member functions
0081 //
0082 
0083 // ------------ method called to produce the data  ------------
0084 void Phase2L1CaloPFClusterEmulator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0085   using namespace edm;
0086   std::unique_ptr<l1tp2::CaloPFClusterCollection> pfclusterCands(make_unique<l1tp2::CaloPFClusterCollection>());
0087 
0088   edm::Handle<std::vector<l1tp2::CaloTower>> caloTowerCollection;
0089   iEvent.getByToken(caloTowerToken_, caloTowerCollection);
0090   if (!caloTowerCollection.isValid())
0091     cms::Exception("Phase2L1CaloPFClusterEmulator") << "Failed to get towers from caloTowerCollection!";
0092 
0093   float GCTintTowers[nTowerEta][nTowerPhi];
0094   float realEta[nTowerEta][nTowerPhi];
0095   float realPhi[nTowerEta][nTowerPhi];
0096   for (const l1tp2::CaloTower& i : *caloTowerCollection) {
0097     int ieta = i.towerIEta();
0098     int iphi = i.towerIPhi();
0099     GCTintTowers[ieta][iphi] = i.ecalTowerEt();
0100     realEta[ieta][iphi] = i.towerEta();
0101     realPhi[ieta][iphi] = i.towerPhi();
0102   }
0103 
0104   float regions[nSLR][nTowerEtaSLR][nTowerPhiSLR];
0105 
0106   for (int ieta = 0; ieta < nTowerEtaSLR; ieta++) {
0107     for (int iphi = 0; iphi < nTowerPhiSLR; iphi++) {
0108       for (int k = 0; k < nSLR; k++) {
0109         regions[k][ieta][iphi] = 0.;
0110       }
0111     }
0112   }
0113 
0114   //Assign ETs to 36 21x8 regions
0115 
0116   for (int ieta = 0; ieta < nTowerEtaSLR; ieta++) {
0117     for (int iphi = 0; iphi < nTowerPhiSLR; iphi++) {
0118       if (ieta > 1) {
0119         if (iphi > 1)
0120           regions[0][ieta][iphi] = GCTintTowers[ieta - 2][iphi - 2];
0121         for (int k = 1; k < 17; k++) {
0122           regions[k * 2][ieta][iphi] = GCTintTowers[ieta - 2][iphi + k * 4 - 2];
0123         }
0124         if (iphi < 6)
0125           regions[34][ieta][iphi] = GCTintTowers[ieta - 2][iphi + 66];
0126       }
0127       if (ieta < 19) {
0128         if (iphi > 1)
0129           regions[1][ieta][iphi] = GCTintTowers[ieta + 15][iphi - 2];
0130         for (int k = 1; k < 17; k++) {
0131           regions[k * 2 + 1][ieta][iphi] = GCTintTowers[ieta + 15][iphi + k * 4 - 2];
0132         }
0133         if (iphi < 6)
0134           regions[35][ieta][iphi] = GCTintTowers[ieta + 15][iphi + 66];
0135       }
0136     }
0137   }
0138 
0139   float temporary[nTowerEtaSLR][nTowerPhiSLR];
0140   int etaoffset = 0;
0141   int phioffset = 0;
0142 
0143   //Use same code from firmware for finding clusters
0144   for (int k = 0; k < nSLR; k++) {
0145     for (int ieta = 0; ieta < nTowerEtaSLR; ieta++) {
0146       for (int iphi = 0; iphi < nTowerPhiSLR; iphi++) {
0147         temporary[ieta][iphi] = regions[k][ieta][iphi];
0148       }
0149     }
0150     if (k % 2 == 0)
0151       etaoffset = 0;
0152     else
0153       etaoffset = 17;
0154     if (k > 1 && k % 2 == 0)
0155       phioffset = phioffset + 4;
0156 
0157     gctpf::PFcluster_t tempPfclusters;
0158     tempPfclusters = gctpf::pfcluster(temporary, etaoffset, phioffset);
0159 
0160     for (int i = 0; i < nPFClusterSLR; i++) {
0161       int gcteta = tempPfclusters.GCTpfclusters[i].eta;
0162       int gctphi = tempPfclusters.GCTpfclusters[i].phi;
0163       float towereta = realEta[gcteta][gctphi];
0164       float towerphi = realPhi[gcteta][gctphi];
0165       l1tp2::CaloPFCluster l1CaloPFCluster;
0166       l1CaloPFCluster.setClusterEt(tempPfclusters.GCTpfclusters[i].et);
0167       l1CaloPFCluster.setClusterIEta(gcteta);
0168       l1CaloPFCluster.setClusterIPhi(gctphi);
0169       l1CaloPFCluster.setClusterEta(towereta);
0170       l1CaloPFCluster.setClusterPhi(towerphi);
0171       pfclusterCands->push_back(l1CaloPFCluster);
0172     }
0173   }
0174 
0175   edm::Handle<HcalTrigPrimDigiCollection> hfHandle;
0176   if (!iEvent.getByToken(hfToken_, hfHandle))
0177     edm::LogError("Phase2L1CaloJetEmulator") << "Failed to get HcalTrigPrimDigi for HF!";
0178   iEvent.getByToken(hfToken_, hfHandle);
0179 
0180   float hfTowers[2 * nHfEta][nHfPhi];  // split 12 -> 24
0181   float hfEta[nHfEta][nHfPhi];
0182   float hfPhi[nHfEta][nHfPhi];
0183   for (int iphi = 0; iphi < nHfPhi; iphi++) {
0184     for (int ieta = 0; ieta < nHfEta; ieta++) {
0185       hfTowers[2 * ieta][iphi] = 0;
0186       hfTowers[2 * ieta + 1][iphi] = 0;
0187       int temp;
0188       if (ieta < nHfEta / 2)
0189         temp = ieta - l1t::CaloTools::kHFEnd;
0190       else
0191         temp = ieta - nHfEta / 2 + l1t::CaloTools::kHFBegin + 1;
0192       hfEta[ieta][iphi] = l1t::CaloTools::towerEta(temp);
0193       hfPhi[ieta][iphi] = l1t::CaloTools::towerPhi(temp, iphi + 1);
0194     }
0195   }
0196 
0197   float regionsHF[nHfRegions][nHfEta]
0198                  [nHfPhi /
0199                   6];  // 6 regions each 3 sectors (1 sector = 20 deg in phi), 12x12->24x12 unique towers, no overlap in phi
0200 
0201   for (int ieta = 0; ieta < nHfEta; ieta++) {
0202     for (int iphi = 0; iphi < nHfPhi / 6; iphi++) {
0203       for (int k = 0; k < nHfRegions; k++) {
0204         regionsHF[k][ieta][iphi] = 0.;
0205       }
0206     }
0207   }
0208 
0209   // Assign ET to 12 24x12 regions
0210   const auto& decoder = iSetup.getData(decoderTag_);
0211   for (const auto& hit : *hfHandle.product()) {
0212     double et = decoder.hcaletValue(hit.id(), hit.t0());
0213     int ieta = 0;
0214     if (abs(hit.id().ieta()) < l1t::CaloTools::kHFBegin)
0215       continue;
0216     if (abs(hit.id().ieta()) > l1t::CaloTools::kHFEnd)
0217       continue;
0218     if (hit.id().ieta() <= -(l1t::CaloTools::kHFBegin + 1)) {
0219       ieta = hit.id().ieta() + l1t::CaloTools::kHFEnd;
0220     } else if (hit.id().ieta() >= (l1t::CaloTools::kHFBegin + 1)) {
0221       ieta = nHfEta / 2 + (hit.id().ieta() - (l1t::CaloTools::kHFBegin + 1));
0222     }
0223     int iphi = hit.id().iphi() - 1;  // HF phi runs between 1-72
0224     // split tower energy
0225     hfTowers[2 * ieta][iphi] = et / 2;
0226     hfTowers[2 * ieta + 1][iphi] = et / 2;
0227     if ((ieta < 2 || ieta >= nHfEta - 2) && iphi % 4 == 2) {
0228       hfTowers[2 * ieta][iphi] = et / 8;
0229       hfTowers[2 * ieta + 1][iphi] = et / 8;
0230       hfTowers[2 * ieta][iphi + 1] = et / 8;
0231       hfTowers[2 * ieta + 1][iphi + 1] = et / 8;
0232       hfTowers[2 * ieta][iphi + 2] = et / 8;
0233       hfTowers[2 * ieta + 1][iphi + 2] = et / 8;
0234       hfTowers[2 * ieta][iphi + 3] = et / 8;
0235       hfTowers[2 * ieta + 1][iphi + 3] = et / 8;
0236     } else if ((ieta >= 2 && ieta < nHfEta - 2) && iphi % 2 == 0) {
0237       hfTowers[2 * ieta][iphi] = et / 4;
0238       hfTowers[2 * ieta + 1][iphi] = et / 4;
0239       hfTowers[2 * ieta][iphi + 1] = et / 4;
0240       hfTowers[2 * ieta + 1][iphi + 1] = et / 4;
0241     }
0242   }
0243 
0244   for (int ieta = 0; ieta < 2 * nHfEta; ieta++) {
0245     for (int iphi = 0; iphi < nHfPhi / 6; iphi++) {
0246       if (ieta < nHfEta) {
0247         regionsHF[0][ieta][0] = hfTowers[ieta][70];
0248         regionsHF[0][ieta][1] = hfTowers[ieta][71];
0249         for (int k = 0; k < nHfRegions; k += 2) {
0250           if (k != 0 || iphi > 1)
0251             regionsHF[k][ieta][iphi] = hfTowers[ieta][iphi + k * (nHfRegions / 2) - 2];
0252         }
0253       } else {
0254         regionsHF[1][ieta - nHfEta][0] = hfTowers[ieta][70];
0255         regionsHF[1][ieta - nHfEta][1] = hfTowers[ieta][71];
0256         for (int k = 1; k < nHfRegions; k += 2) {
0257           if (k != 1 || iphi > 1)
0258             regionsHF[k][ieta - nHfEta][iphi] = hfTowers[ieta][iphi + (k - 1) * (nHfRegions / 2) - 2];
0259         }
0260       }
0261     }
0262   }
0263 
0264   float temporaryHF[nHfEta][nHfPhi / 6];
0265   int etaoffsetHF = 0;
0266   int phioffsetHF = -2;
0267 
0268   for (int k = 0; k < nHfRegions; k++) {
0269     for (int ieta = 0; ieta < nHfEta; ieta++) {
0270       for (int iphi = 0; iphi < nHfPhi / 6; iphi++) {
0271         temporaryHF[ieta][iphi] = regionsHF[k][ieta][iphi];
0272       }
0273     }
0274     if (k % 2 == 0)
0275       etaoffsetHF = 0;
0276     else
0277       etaoffsetHF = nHfEta;
0278     if (k > 1 && k % 2 == 0)
0279       phioffsetHF = phioffsetHF + nHfPhi / 6;
0280 
0281     gctpf::PFcluster_t tempPfclustersHF;
0282     tempPfclustersHF = gctpf::pfclusterHF(temporaryHF, etaoffsetHF, phioffsetHF);
0283 
0284     for (int i = 0; i < nPFClusterSLR; i++) {
0285       int hfeta = tempPfclustersHF.GCTpfclusters[i].eta / 2;  // turn back 24 -> 12
0286       int hfphi = tempPfclustersHF.GCTpfclusters[i].phi;
0287       float towereta = hfEta[hfeta][hfphi];
0288       float towerphi = hfPhi[hfeta][hfphi];
0289       l1tp2::CaloPFCluster l1CaloPFCluster;
0290       l1CaloPFCluster.setClusterEt(tempPfclustersHF.GCTpfclusters[i].et);
0291       l1CaloPFCluster.setClusterIEta(hfeta);
0292       l1CaloPFCluster.setClusterIPhi(hfphi);
0293       l1CaloPFCluster.setClusterEta(towereta);
0294       l1CaloPFCluster.setClusterPhi(towerphi);
0295       pfclusterCands->push_back(l1CaloPFCluster);
0296     }
0297   }
0298 
0299   iEvent.put(std::move(pfclusterCands), "GCTPFCluster");
0300 }
0301 
0302 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0303 void Phase2L1CaloPFClusterEmulator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0304   edm::ParameterSetDescription desc;
0305   desc.add<edm::InputTag>("gctFullTowers", edm::InputTag("l1tPhase2L1CaloEGammaEmulator", "GCTFullTowers"));
0306   desc.add<edm::InputTag>("hcalDigis", edm::InputTag("simHcalTriggerPrimitiveDigis"));
0307   descriptions.addWithDefaultLabel(desc);
0308 }
0309 
0310 //define this as a plug-in
0311 DEFINE_FWK_MODULE(Phase2L1CaloPFClusterEmulator);