File indexing completed on 2024-09-05 05:04:37
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021
0022
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
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
0064 edm::EDGetTokenT<l1tp2::CaloTowerCollection> caloTowerToken_;
0065 edm::EDGetTokenT<HcalTrigPrimDigiCollection> hfToken_;
0066 edm::ESGetToken<CaloTPGTranscoder, CaloTPGRecord> decoderTag_;
0067 };
0068
0069
0070
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
0081
0082
0083
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
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
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];
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];
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
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;
0224
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;
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
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
0311 DEFINE_FWK_MODULE(Phase2L1CaloPFClusterEmulator);