Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-12-20 03:13:54

0001 // -*- C++ -*-
0002 //
0003 // Package:    L1Trigger/L1TCaloSummary
0004 // Class:      L1TCaloSummary
0005 //
0006 /**\class L1TCaloSummary L1TCaloSummary.cc L1Trigger/L1TCaloSummary/plugins/L1TCaloSummary.cc
0007 
0008    Description: The package L1Trigger/L1TCaloSummary is prepared for monitoring the CMS Layer-1 Calorimeter Trigger.
0009 
0010    Implementation:
0011    It prepares region objects and puts them in the event
0012 */
0013 //
0014 // Original Author:  Sridhara Dasu
0015 //         Created:  Sat, 14 Nov 2015 14:18:27 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 
0031 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0032 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0033 
0034 #include "L1Trigger/L1TCaloLayer1/src/UCTLayer1.hh"
0035 #include "L1Trigger/L1TCaloLayer1/src/UCTCrate.hh"
0036 #include "L1Trigger/L1TCaloLayer1/src/UCTCard.hh"
0037 #include "L1Trigger/L1TCaloLayer1/src/UCTRegion.hh"
0038 #include "L1Trigger/L1TCaloLayer1/src/UCTTower.hh"
0039 #include "L1Trigger/L1TCaloLayer1/src/UCTGeometry.hh"
0040 
0041 #include "L1Trigger/L1TCaloLayer1/src/UCTObject.hh"
0042 #include "L1Trigger/L1TCaloLayer1/src/UCTSummaryCard.hh"
0043 #include "L1Trigger/L1TCaloLayer1/src/UCTGeometryExtended.hh"
0044 
0045 #include "DataFormats/L1Trigger/interface/L1EmParticle.h"
0046 #include "DataFormats/L1Trigger/interface/L1EmParticleFwd.h"
0047 #include "DataFormats/L1Trigger/interface/L1JetParticle.h"
0048 #include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h"
0049 #include "DataFormats/L1Trigger/interface/L1EtMissParticle.h"
0050 #include "DataFormats/L1Trigger/interface/L1EtMissParticleFwd.h"
0051 
0052 #include "DataFormats/L1CaloTrigger/interface/L1CaloCollections.h"
0053 #include "DataFormats/L1CaloTrigger/interface/L1CaloRegion.h"
0054 #include "DataFormats/L1CaloTrigger/interface/CICADA.h"
0055 
0056 #include "DataFormats/Math/interface/LorentzVector.h"
0057 
0058 #include "L1Trigger/L1TCaloLayer1/src/UCTLogging.hh"
0059 #include <bitset>
0060 
0061 #include <string>
0062 #include <sstream>
0063 
0064 //Anomaly detection includes
0065 #include "ap_fixed.h"
0066 #include "hls4ml/emulator.h"
0067 
0068 using namespace l1tcalo;
0069 using namespace l1extra;
0070 using namespace std;
0071 
0072 //
0073 // class declaration
0074 //
0075 
0076 template <class INPUT, class OUTPUT>
0077 class L1TCaloSummary : public edm::stream::EDProducer<> {
0078 public:
0079   explicit L1TCaloSummary(const edm::ParameterSet&);
0080   ~L1TCaloSummary() override = default;
0081 
0082   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0083 
0084 private:
0085   //void beginJob() override;
0086   void produce(edm::Event&, const edm::EventSetup&) override;
0087   //void endJob() override;
0088 
0089   void print();
0090 
0091   // ----------member data ---------------------------
0092 
0093   uint32_t nPumBins;
0094 
0095   std::vector<std::vector<std::vector<uint32_t>>> pumLUT;
0096 
0097   double caloScaleFactor;
0098 
0099   uint32_t jetSeed;
0100   uint32_t tauSeed;
0101   float tauIsolationFactor;
0102   uint32_t eGammaSeed;
0103   double eGammaIsolationFactor;
0104   double boostedJetPtFactor;
0105 
0106   bool verbose;
0107   int fwVersion;
0108 
0109   edm::EDGetTokenT<L1CaloRegionCollection> regionToken;
0110   edm::EDGetTokenT<L1CaloRegionCollection> backupRegionToken;
0111 
0112   UCTLayer1* layer1;
0113 
0114   hls4mlEmulator::ModelLoader loader;
0115   std::shared_ptr<hls4mlEmulator::Model> model;
0116 
0117   bool overwriteWithTestPatterns;
0118   std::vector<edm::ParameterSet> testPatterns;
0119 };
0120 
0121 //
0122 // constants, enums and typedefs
0123 //
0124 
0125 //
0126 // static data member definitions
0127 //
0128 
0129 //
0130 // constructors and destructor
0131 //
0132 template <class INPUT, class OUTPUT>
0133 L1TCaloSummary<INPUT, OUTPUT>::L1TCaloSummary(const edm::ParameterSet& iConfig)
0134     : nPumBins(iConfig.getParameter<unsigned int>("nPumBins")),
0135       pumLUT(nPumBins, std::vector<std::vector<uint32_t>>(2, std::vector<uint32_t>(13))),
0136       caloScaleFactor(iConfig.getParameter<double>("caloScaleFactor")),
0137       jetSeed(iConfig.getParameter<unsigned int>("jetSeed")),
0138       tauSeed(iConfig.getParameter<unsigned int>("tauSeed")),
0139       tauIsolationFactor(iConfig.getParameter<double>("tauIsolationFactor")),
0140       eGammaSeed(iConfig.getParameter<unsigned int>("eGammaSeed")),
0141       eGammaIsolationFactor(iConfig.getParameter<double>("eGammaIsolationFactor")),
0142       boostedJetPtFactor(iConfig.getParameter<double>("boostedJetPtFactor")),
0143       verbose(iConfig.getParameter<bool>("verbose")),
0144       fwVersion(iConfig.getParameter<int>("firmwareVersion")),
0145       regionToken(consumes<L1CaloRegionCollection>(iConfig.getParameter<edm::InputTag>("caloLayer1Regions"))),
0146       //backupRegionToken(consumes<L1CaloRegionCollection>(edm::InputTag("simCaloStage2Layer1Digis"))),
0147       backupRegionToken(consumes<L1CaloRegionCollection>(iConfig.getParameter<edm::InputTag>("backupRegionToken"))),
0148       loader(hls4mlEmulator::ModelLoader(iConfig.getParameter<string>("CICADAModelVersion"))),
0149       overwriteWithTestPatterns(iConfig.getParameter<bool>("useTestPatterns")),
0150       testPatterns(iConfig.getParameter<std::vector<edm::ParameterSet>>("testPatterns")) {
0151   std::vector<double> pumLUTData;
0152   char pumLUTString[10];
0153   for (uint32_t pumBin = 0; pumBin < nPumBins; pumBin++) {
0154     for (uint32_t side = 0; side < 2; side++) {
0155       if (side == 0)
0156         sprintf(pumLUTString, "pumLUT%2.2dp", pumBin);
0157       else
0158         sprintf(pumLUTString, "pumLUT%2.2dn", pumBin);
0159       pumLUTData = iConfig.getParameter<std::vector<double>>(pumLUTString);
0160       for (uint32_t iEta = 0; iEta < std::max((uint32_t)pumLUTData.size(), MaxUCTRegionsEta); iEta++) {
0161         pumLUT[pumBin][side][iEta] = (uint32_t)round(pumLUTData[iEta] / caloScaleFactor);
0162       }
0163       if (pumLUTData.size() != (MaxUCTRegionsEta))
0164         edm::LogError("L1TCaloSummary") << "PUM LUT Data size integrity check failed; Expected size = "
0165                                         << MaxUCTRegionsEta << "; Provided size = " << pumLUTData.size()
0166                                         << "; Will use what is provided :(" << std::endl;
0167     }
0168   }
0169   produces<L1JetParticleCollection>("Boosted");
0170 
0171   //anomaly trigger loading
0172   model = loader.load_model();
0173   produces<l1t::CICADABxCollection>("CICADAScore");
0174 }
0175 
0176 //
0177 // member functions
0178 //
0179 
0180 // ------------ method called to produce the data  ------------
0181 template <class INPUT, class OUTPUT>
0182 void L1TCaloSummary<INPUT, OUTPUT>::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0183   using namespace edm;
0184 
0185   std::unique_ptr<L1JetParticleCollection> bJetCands(new L1JetParticleCollection);
0186 
0187   std::unique_ptr<l1t::CICADABxCollection> CICADAScore = std::make_unique<l1t::CICADABxCollection>();
0188   CICADAScore->setBXRange(-2, 2);
0189 
0190   UCTGeometry g;
0191 
0192   // Here we read region data from the region collection created by L1TCaloLayer1 instead of
0193   // independently creating regions from TPGs for processing by the summary card. This results
0194   // in a single region vector of size 252 whereas from independent creation we had 3*6 vectors
0195   // of size 7*2. Indices are mapped in UCTSummaryCard accordingly.
0196   UCTSummaryCard summaryCard =
0197       UCTSummaryCard(&pumLUT, jetSeed, tauSeed, tauIsolationFactor, eGammaSeed, eGammaIsolationFactor);
0198   std::vector<UCTRegion> inputRegions;
0199   inputRegions.reserve(252);  // 252 calorimeter regions. 18 phi * 14 eta
0200   edm::Handle<std::vector<L1CaloRegion>> regionCollection;
0201   if (!iEvent.getByToken(regionToken, regionCollection))
0202     edm::LogError("L1TCaloSummary") << "UCT: Failed to get regions from region collection!";
0203   iEvent.getByToken(regionToken, regionCollection);
0204 
0205   if (regionCollection->empty()) {
0206     iEvent.getByToken(backupRegionToken, regionCollection);
0207     edm::LogWarning("L1TCaloSummary") << "Switched to emulated regions since data regions was empty.\n";
0208   }
0209 
0210   //Model input
0211   //This is done as a flat vector input, but future versions may involve 2D input
0212   //This will have to be handled later
0213   INPUT modelInput[252];
0214   for (const L1CaloRegion& i : *regionCollection) {
0215     UCTRegionIndex r = g.getUCTRegionIndexFromL1CaloRegion(i.gctEta(), i.gctPhi());
0216     UCTTowerIndex t = g.getUCTTowerIndexFromL1CaloRegion(r, i.raw());
0217     uint32_t absCaloEta = std::abs(t.first);
0218     uint32_t absCaloPhi = std::abs(t.second);
0219     bool negativeEta = false;
0220     if (t.first < 0)
0221       negativeEta = true;
0222     uint32_t crate = g.getCrate(t.first, t.second);
0223     uint32_t card = g.getCard(t.first, t.second);
0224     uint32_t region = g.getRegion(absCaloEta, absCaloPhi);
0225     UCTRegion test = UCTRegion(crate, card, negativeEta, region, fwVersion);
0226     test.setRegionSummary(i.raw());
0227     inputRegions.push_back(test);
0228     //This *should* fill the tensor in the proper order to be fed to the anomaly model
0229     //We take 4 off of the GCT eta/iEta.
0230     //iEta taken from this ranges from 4-17, (I assume reserving lower and higher for forward regions)
0231     //So our first index, index 0, is technically iEta=4, and so-on.
0232     //CICADA reads this as a flat vector
0233     modelInput[14 * i.gctPhi() + (i.gctEta() - 4)] = i.et();
0234   }
0235   // Check if we're using test patterns. If so, we overwrite the inputs with a test pattern
0236   if (overwriteWithTestPatterns) {
0237     unsigned int evt = iEvent.id().event();
0238     unsigned int totalTestPatterns = testPatterns.size();
0239     unsigned int patternElement = evt % totalTestPatterns;
0240     const edm::ParameterSet& element = testPatterns.at(patternElement);
0241     std::stringstream inputStream;
0242     std::string PhiRowString;
0243 
0244     edm::LogWarning("L1TCaloSummary") << "Overwriting existing CICADA input with test pattern!\n";
0245 
0246     for (unsigned short int iPhi = 1; iPhi <= 18; ++iPhi) {
0247       PhiRowString = "";
0248       std::stringstream PhiRowStringStream;
0249       PhiRowStringStream << "iPhi_" << iPhi;
0250       PhiRowString = PhiRowStringStream.str();
0251       std::vector<unsigned int> phiRow = element.getParameter<std::vector<unsigned int>>(PhiRowString);
0252       for (unsigned short int iEta = 1; iEta <= 14; ++iEta) {
0253         modelInput[14 * (iPhi - 1) + (iEta - 1)] = phiRow.at(iEta - 1);
0254         inputStream << phiRow.at(iEta - 1) << " ";
0255       }
0256       inputStream << "\n";
0257     }
0258     edm::LogInfo("L1TCaloSummary") << "Input Stream:\n" << inputStream.str();
0259   }
0260 
0261   //Extract model output
0262   OUTPUT modelResult[1] = {
0263       OUTPUT("0.0", 10)};  //the 10 here refers to the fact that we read in "0.0" as a decimal number
0264   model->prepare_input(modelInput);
0265   model->predict();
0266   model->read_result(modelResult);
0267 
0268   CICADAScore->push_back(0, modelResult[0].to_float());
0269 
0270   if (overwriteWithTestPatterns)
0271     edm::LogInfo("L1TCaloSummary") << "Test Pattern Output: " << CICADAScore->at(0, 0);
0272 
0273   summaryCard.setRegionData(inputRegions);
0274 
0275   if (!summaryCard.process()) {
0276     edm::LogError("L1TCaloSummary") << "UCT: Failed to process summary card" << std::endl;
0277     exit(1);
0278   }
0279 
0280   double pt = 0;
0281   double eta = -999.;
0282   double phi = -999.;
0283   double mass = 0;
0284 
0285   std::list<std::shared_ptr<UCTObject>> boostedJetObjs = summaryCard.getBoostedJetObjs();
0286   for (std::list<std::shared_ptr<UCTObject>>::const_iterator i = boostedJetObjs.begin(); i != boostedJetObjs.end();
0287        i++) {
0288     const std::shared_ptr<UCTObject>& object = *i;
0289     pt = ((double)object->et()) * caloScaleFactor * boostedJetPtFactor;
0290     eta = g.getUCTTowerEta(object->iEta());
0291     phi = g.getUCTTowerPhi(object->iPhi());
0292     bitset<3> activeRegionEtaPattern = 0;
0293     for (uint32_t iEta = 0; iEta < 3; iEta++) {
0294       bool activeStrip = false;
0295       for (uint32_t iPhi = 0; iPhi < 3; iPhi++) {
0296         if (object->boostedJetRegionET()[3 * iEta + iPhi] > 30 &&
0297             object->boostedJetRegionET()[3 * iEta + iPhi] > object->et() * 0.0625)
0298           activeStrip = true;
0299       }
0300       if (activeStrip)
0301         activeRegionEtaPattern |= (0x1 << iEta);
0302     }
0303     bitset<3> activeRegionPhiPattern = 0;
0304     for (uint32_t iPhi = 0; iPhi < 3; iPhi++) {
0305       bool activeStrip = false;
0306       for (uint32_t iEta = 0; iEta < 3; iEta++) {
0307         if (object->boostedJetRegionET()[3 * iEta + iPhi] > 30 &&
0308             object->boostedJetRegionET()[3 * iEta + iPhi] > object->et() * 0.0625)
0309           activeStrip = true;
0310       }
0311       if (activeStrip)
0312         activeRegionPhiPattern |= (0x1 << iPhi);
0313     }
0314     string regionEta = activeRegionEtaPattern.to_string<char, std::string::traits_type, std::string::allocator_type>();
0315     string regionPhi = activeRegionPhiPattern.to_string<char, std::string::traits_type, std::string::allocator_type>();
0316 
0317     bool centralHighest = object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[0] &&
0318                           object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[1] &&
0319                           object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[2] &&
0320                           object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[3] &&
0321                           object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[5] &&
0322                           object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[6] &&
0323                           object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[7] &&
0324                           object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[8];
0325 
0326     if (abs(eta) < 2.5 && ((regionEta == "101" && (regionPhi == "110" || regionPhi == "101" || regionPhi == "010")) ||
0327                            ((regionEta == "110" || regionEta == "101" || regionEta == "010") && regionPhi == "101") ||
0328                            (regionEta == "111" && (regionPhi == "110" || regionPhi == "010")) ||
0329                            ((regionEta == "110" || regionEta == "010") && regionPhi == "111") ||
0330                            ((regionEta == "010" || regionPhi == "010" || regionEta == "110" || regionPhi == "110" ||
0331                              regionEta == "011" || regionPhi == "011") &&
0332                             centralHighest)))
0333       bJetCands->push_back(L1JetParticle(math::PtEtaPhiMLorentzVector(pt, eta, phi, mass), L1JetParticle::kCentral));
0334   }
0335 
0336   iEvent.put(std::move(bJetCands), "Boosted");
0337   //Write out anomaly score
0338   iEvent.put(std::move(CICADAScore), "CICADAScore");
0339 }
0340 
0341 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0342 template <class INPUT, class OUTPUT>
0343 void L1TCaloSummary<INPUT, OUTPUT>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0344   //The following says we do not know what parameters are allowed so do no validation
0345   // Please change this to state exactly what you do use, even if it is no parameters
0346   edm::ParameterSetDescription desc;
0347   desc.setUnknown();
0348   descriptions.addDefault(desc);
0349 }
0350 
0351 // Initial version, X.0.0, input/output typing
0352 typedef L1TCaloSummary<ap_ufixed<10, 10>, ap_fixed<11, 5>> L1TCaloSummary_CICADA_v1p0p0;
0353 typedef L1TCaloSummary<ap_uint<10>, ap_ufixed<16, 8>> L1TCaloSummary_CICADA_v2p0p0;
0354 DEFINE_FWK_MODULE(L1TCaloSummary_CICADA_v1p0p0);
0355 DEFINE_FWK_MODULE(L1TCaloSummary_CICADA_v2p0p0);
0356 // X.1.0 version input.output typing
0357 typedef L1TCaloSummary<ap_ufixed<10, 10>, ap_fixed<11, 5>> L1TCaloSummary_CICADA_v1p1p0;
0358 typedef L1TCaloSummary<ap_uint<10>, ap_ufixed<16, 8>> L1TCaloSummary_CICADA_v2p1p0;
0359 DEFINE_FWK_MODULE(L1TCaloSummary_CICADA_v1p1p0);
0360 DEFINE_FWK_MODULE(L1TCaloSummary_CICADA_v2p1p0);
0361 // X.1.1 version input/output typing
0362 typedef L1TCaloSummary<ap_uint<10>, ap_ufixed<16, 8, AP_RND, AP_SAT, AP_SAT>> L1TCaloSummary_CICADA_vXp1p1;
0363 DEFINE_FWK_MODULE(L1TCaloSummary_CICADA_vXp1p1);
0364 // X.1.2 version input/output typing
0365 typedef L1TCaloSummary<ap_uint<10>, ap_ufixed<16, 8, AP_RND_CONV, AP_SAT>> L1TCaloSummary_CICADA_vXp1p2;
0366 DEFINE_FWK_MODULE(L1TCaloSummary_CICADA_vXp1p2);