L1TCaloLayer1/plugins/L1TCaloSummary.cc

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
0055 #include "DataFormats/Math/interface/LorentzVector.h"
0056
0057 #include "L1Trigger/L1TCaloLayer1/src/UCTLogging.hh"
0058 #include <bitset>
0059
0060 //Anomaly detection includes
0061 #include "ap_fixed.h"
0062 #include "hls4ml/emulator.h"
0063
0064 using namespace l1tcalo;
0065 using namespace l1extra;
0066 using namespace std;
0067
0068 //
0069 // class declaration
0070 //
0071
0072 template <class INPUT, class OUTPUT>
0073 class L1TCaloSummary : public edm::stream::EDProducer<> {
0074 public:
0075   explicit L1TCaloSummary(const edm::ParameterSet&);
0076   ~L1TCaloSummary() override = default;
0077
0078   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0079
0080 private:
0081   //void beginJob() override;
0082   void produce(edm::Event&, const edm::EventSetup&) override;
0083   //void endJob() override;
0084
0085   void beginRun(edm::Run const&, edm::EventSetup const&) override{};
0086
0087   void print();
0088
0089   // ----------member data ---------------------------
0090
0091   uint32_t nPumBins;
0092
0093   std::vector<std::vector<std::vector<uint32_t>>> pumLUT;
0094
0095   double caloScaleFactor;
0096
0097   uint32_t jetSeed;
0098   uint32_t tauSeed;
0099   float tauIsolationFactor;
0100   uint32_t eGammaSeed;
0101   double eGammaIsolationFactor;
0102   double boostedJetPtFactor;
0103
0104   bool verbose;
0105   int fwVersion;
0106
0107   edm::EDGetTokenT<L1CaloRegionCollection> regionToken;
0108
0109   UCTLayer1* layer1;
0110
0111   hls4mlEmulator::ModelLoader loader;
0112   std::shared_ptr<hls4mlEmulator::Model> model;
0113 };
0114
0115 //
0116 // constants, enums and typedefs
0117 //
0118
0119 //
0120 // static data member definitions
0121 //
0122
0123 //
0124 // constructors and destructor
0125 //
0126 template <class INPUT, class OUTPUT>
0127 L1TCaloSummary<INPUT, OUTPUT>::L1TCaloSummary(const edm::ParameterSet& iConfig)
0128     : nPumBins(iConfig.getParameter<unsigned int>("nPumBins")),
0129       pumLUT(nPumBins, std::vector<std::vector<uint32_t>>(2, std::vector<uint32_t>(13))),
0130       caloScaleFactor(iConfig.getParameter<double>("caloScaleFactor")),
0131       jetSeed(iConfig.getParameter<unsigned int>("jetSeed")),
0132       tauSeed(iConfig.getParameter<unsigned int>("tauSeed")),
0133       tauIsolationFactor(iConfig.getParameter<double>("tauIsolationFactor")),
0134       eGammaSeed(iConfig.getParameter<unsigned int>("eGammaSeed")),
0135       eGammaIsolationFactor(iConfig.getParameter<double>("eGammaIsolationFactor")),
0136       boostedJetPtFactor(iConfig.getParameter<double>("boostedJetPtFactor")),
0137       verbose(iConfig.getParameter<bool>("verbose")),
0138       fwVersion(iConfig.getParameter<int>("firmwareVersion")),
0139       regionToken(consumes<L1CaloRegionCollection>(edm::InputTag("simCaloStage2Layer1Digis"))),
0140       loader(hls4mlEmulator::ModelLoader(iConfig.getParameter<string>("CICADAModelVersion"))) {
0141   std::vector<double> pumLUTData;
0142   char pumLUTString[10];
0143   for (uint32_t pumBin = 0; pumBin < nPumBins; pumBin++) {
0144     for (uint32_t side = 0; side < 2; side++) {
0145       if (side == 0)
0146         sprintf(pumLUTString, "pumLUT%2.2dp", pumBin);
0147       else
0148         sprintf(pumLUTString, "pumLUT%2.2dn", pumBin);
0149       pumLUTData = iConfig.getParameter<std::vector<double>>(pumLUTString);
0150       for (uint32_t iEta = 0; iEta < std::max((uint32_t)pumLUTData.size(), MaxUCTRegionsEta); iEta++) {
0151         pumLUT[pumBin][side][iEta] = (uint32_t)round(pumLUTData[iEta] / caloScaleFactor);
0152       }
0153       if (pumLUTData.size() != (MaxUCTRegionsEta))
0154         edm::LogError("L1TCaloSummary") << "PUM LUT Data size integrity check failed; Expected size = "
0155                                         << MaxUCTRegionsEta << "; Provided size = " << pumLUTData.size()
0156                                         << "; Will use what is provided :(" << std::endl;
0157     }
0158   }
0159   produces<L1JetParticleCollection>("Boosted");
0160
0161   //anomaly trigger loading
0162   model = loader.load_model();
0163   produces<float>("CICADAScore");
0164 }
0165
0166 //
0167 // member functions
0168 //
0169
0170 // ------------ method called to produce the data  ------------
0171 template <class INPUT, class OUTPUT>
0172 void L1TCaloSummary<INPUT, OUTPUT>::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0173   using namespace edm;
0174
0175   std::unique_ptr<L1JetParticleCollection> bJetCands(new L1JetParticleCollection);
0176
0177   std::unique_ptr<float> CICADAScore = std::make_unique<float>();
0178
0179   UCTGeometry g;
0180
0181   // Here we read region data from the region collection created by L1TCaloLayer1 instead of
0182   // independently creating regions from TPGs for processing by the summary card. This results
0183   // in a single region vector of size 252 whereas from independent creation we had 3*6 vectors
0184   // of size 7*2. Indices are mapped in UCTSummaryCard accordingly.
0185   UCTSummaryCard summaryCard =
0186       UCTSummaryCard(&pumLUT, jetSeed, tauSeed, tauIsolationFactor, eGammaSeed, eGammaIsolationFactor);
0187   std::vector<UCTRegion*> inputRegions;
0188   inputRegions.clear();
0189   edm::Handle<std::vector<L1CaloRegion>> regionCollection;
0190   if (!iEvent.getByToken(regionToken, regionCollection))
0191     edm::LogError("L1TCaloSummary") << "UCT: Failed to get regions from region collection!";
0192   iEvent.getByToken(regionToken, regionCollection);
0193   //Model input
0194   //This is done as a flat vector input, but future versions may involve 2D input
0195   //This will have to be handled later
0196   INPUT modelInput[252];
0197   for (const L1CaloRegion& i : *regionCollection) {
0198     UCTRegionIndex r = g.getUCTRegionIndexFromL1CaloRegion(i.gctEta(), i.gctPhi());
0199     UCTTowerIndex t = g.getUCTTowerIndexFromL1CaloRegion(r, i.raw());
0200     uint32_t absCaloEta = std::abs(t.first);
0201     uint32_t absCaloPhi = std::abs(t.second);
0202     bool negativeEta = false;
0203     if (t.first < 0)
0204       negativeEta = true;
0205     uint32_t crate = g.getCrate(t.first, t.second);
0206     uint32_t card = g.getCard(t.first, t.second);
0207     uint32_t region = g.getRegion(absCaloEta, absCaloPhi);
0208     UCTRegion* test = new UCTRegion(crate, card, negativeEta, region, fwVersion);
0209     test->setRegionSummary(i.raw());
0210     inputRegions.push_back(test);
0211     //This *should* fill the tensor in the proper order to be fed to the anomaly model
0212     //We take 4 off of the GCT eta/iEta.
0213     //iEta taken from this ranges from 4-17, (I assume reserving lower and higher for forward regions)
0214     //So our first index, index 0, is technically iEta=4, and so-on.
0215     //CICADA reads this as a flat vector
0216     modelInput[14 * i.gctPhi() + (i.gctEta() - 4)] = i.et();
0217   }
0218   //Extract model output
0219   OUTPUT modelResult[1] = {
0220       OUTPUT("0.0", 10)};  //the 10 here refers to the fact that we read in "0.0" as a decimal number
0221   model->prepare_input(modelInput);
0222   model->predict();
0223   model->read_result(modelResult);
0224
0225   *CICADAScore = modelResult[0].to_float();
0226
0227   summaryCard.setRegionData(inputRegions);
0228
0229   if (!summaryCard.process()) {
0230     edm::LogError("L1TCaloSummary") << "UCT: Failed to process summary card" << std::endl;
0231     exit(1);
0232   }
0233
0234   double pt = 0;
0235   double eta = -999.;
0236   double phi = -999.;
0237   double mass = 0;
0238
0239   std::list<UCTObject*> boostedJetObjs = summaryCard.getBoostedJetObjs();
0240   for (std::list<UCTObject*>::const_iterator i = boostedJetObjs.begin(); i != boostedJetObjs.end(); i++) {
0241     const UCTObject* object = *i;
0242     pt = ((double)object->et()) * caloScaleFactor * boostedJetPtFactor;
0243     eta = g.getUCTTowerEta(object->iEta());
0244     phi = g.getUCTTowerPhi(object->iPhi());
0245     bitset<3> activeRegionEtaPattern = 0;
0246     for (uint32_t iEta = 0; iEta < 3; iEta++) {
0247       bool activeStrip = false;
0248       for (uint32_t iPhi = 0; iPhi < 3; iPhi++) {
0249         if (object->boostedJetRegionET()[3 * iEta + iPhi] > 30 &&
0250             object->boostedJetRegionET()[3 * iEta + iPhi] > object->et() * 0.0625)
0251           activeStrip = true;
0252       }
0253       if (activeStrip)
0254         activeRegionEtaPattern |= (0x1 << iEta);
0255     }
0256     bitset<3> activeRegionPhiPattern = 0;
0257     for (uint32_t iPhi = 0; iPhi < 3; iPhi++) {
0258       bool activeStrip = false;
0259       for (uint32_t iEta = 0; iEta < 3; iEta++) {
0260         if (object->boostedJetRegionET()[3 * iEta + iPhi] > 30 &&
0261             object->boostedJetRegionET()[3 * iEta + iPhi] > object->et() * 0.0625)
0262           activeStrip = true;
0263       }
0264       if (activeStrip)
0265         activeRegionPhiPattern |= (0x1 << iPhi);
0266     }
0267     string regionEta = activeRegionEtaPattern.to_string<char, std::string::traits_type, std::string::allocator_type>();
0268     string regionPhi = activeRegionPhiPattern.to_string<char, std::string::traits_type, std::string::allocator_type>();
0269
0270     bool centralHighest = object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[0] &&
0271                           object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[1] &&
0272                           object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[2] &&
0273                           object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[3] &&
0274                           object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[5] &&
0275                           object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[6] &&
0276                           object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[7] &&
0277                           object->boostedJetRegionET()[4] >= object->boostedJetRegionET()[8];
0278
0279     if (abs(eta) < 2.5 && ((regionEta == "101" && (regionPhi == "110" || regionPhi == "101" || regionPhi == "010")) ||
0280                            ((regionEta == "110" || regionEta == "101" || regionEta == "010") && regionPhi == "101") ||
0281                            (regionEta == "111" && (regionPhi == "110" || regionPhi == "010")) ||
0282                            ((regionEta == "110" || regionEta == "010") && regionPhi == "111") ||
0283                            ((regionEta == "010" || regionPhi == "010" || regionEta == "110" || regionPhi == "110" ||
0284                              regionEta == "011" || regionPhi == "011") &&
0285                             centralHighest)))
0286       bJetCands->push_back(L1JetParticle(math::PtEtaPhiMLorentzVector(pt, eta, phi, mass), L1JetParticle::kCentral));
0287   }
0288
0289   iEvent.put(std::move(bJetCands), "Boosted");
0290   //Write out anomaly score
0291   iEvent.put(std::move(CICADAScore), "CICADAScore");
0292 }
0293
0294 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0295 template <class INPUT, class OUTPUT>
0296 void L1TCaloSummary<INPUT, OUTPUT>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0297   //The following says we do not know what parameters are allowed so do no validation
0298   // Please change this to state exactly what you do use, even if it is no parameters
0299   edm::ParameterSetDescription desc;
0300   desc.setUnknown();
0301   descriptions.addDefault(desc);
0302 }
0303
0304 typedef L1TCaloSummary<ap_ufixed<10, 10>, ap_fixed<11, 5>> L1TCaloSummaryCICADAv1;
0305 typedef L1TCaloSummary<ap_uint<10>, ap_ufixed<16, 8>> L1TCaloSummaryCICADAv2;
0306 //define type version plugins
0307 DEFINE_FWK_MODULE(L1TCaloSummaryCICADAv1);
0308 DEFINE_FWK_MODULE(L1TCaloSummaryCICADAv2);