From 4f0ddb19a673bf3875956f454e4eb05b2b78ce24 Mon Sep 17 00:00:00 2001 From: Peter Eduard Meiring Date: Wed, 29 Jun 2022 10:13:55 +0200 Subject: [PATCH 01/13] Modify emulator scheleton for integration of Composite ID. Adapt the data-formats and preselect candidates for BDT evaluation. Composite ID can be run as alternative to Elliptic Matching. --- .../L1TCorrelator/interface/TkElectron.h | 5 +- DataFormats/L1TCorrelator/src/classes_def.xml | 3 +- .../L1TParticleFlow/interface/PFCandidate.h | 7 +- .../interface/layer1_emulator.h | 2 + .../L1TParticleFlow/src/PFCandidate.cc | 2 + .../L1TParticleFlow/src/classes_def.xml | 4 +- .../HGCalTriggerNtupleHGCMulticlusters.cc | 12 +- .../interface/egamma/pftkegalgo_ref.h | 87 ++++++++- .../plugins/L1TCorrelatorLayer1Producer.cc | 45 +++++ .../python/l1TkEgAlgoEmulator_cfi.py | 27 ++- .../python/l1ctLayer1_cff.py | 3 +- .../src/egamma/pftkegalgo_ref.cpp | 181 +++++++++++++++++- .../src/newfirmware/egamma/compositeID.json | 1 + .../src/newfirmware/egamma/conifer.h | 158 +++++++++++++++ .../test/make_l1ctLayer1_dumpFiles_cfg.py | 2 +- 15 files changed, 509 insertions(+), 30 deletions(-) create mode 100644 L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/compositeID.json create mode 100644 L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/conifer.h diff --git a/DataFormats/L1TCorrelator/interface/TkElectron.h b/DataFormats/L1TCorrelator/interface/TkElectron.h index 9f1ac1eaeeff0..9d66291cbc939 100644 --- a/DataFormats/L1TCorrelator/interface/TkElectron.h +++ b/DataFormats/L1TCorrelator/interface/TkElectron.h @@ -39,16 +39,19 @@ namespace l1t { float trkzVtx() const { return trkzVtx_; } double trackCurvature() const { return trackCurvature_; } - + float compositeBdtScore() const { return compositeBdtScore_; } // ---------- member functions --------------------------- void setTrkzVtx(float TrkzVtx) { trkzVtx_ = TrkzVtx; } void setTrackCurvature(double trackCurvature) { trackCurvature_ = trackCurvature; } + void setCompositeBdtScore(float score) { compositeBdtScore_ = score; } + private: edm::Ptr trkPtr_; float trkzVtx_; double trackCurvature_; + float compositeBdtScore_; }; } // namespace l1t #endif diff --git a/DataFormats/L1TCorrelator/src/classes_def.xml b/DataFormats/L1TCorrelator/src/classes_def.xml index 0685c6bc1356e..47267dc2b8eaa 100644 --- a/DataFormats/L1TCorrelator/src/classes_def.xml +++ b/DataFormats/L1TCorrelator/src/classes_def.xml @@ -44,7 +44,8 @@ - + + diff --git a/DataFormats/L1TParticleFlow/interface/PFCandidate.h b/DataFormats/L1TParticleFlow/interface/PFCandidate.h index c6d52246a6fa5..50cf529178dc9 100644 --- a/DataFormats/L1TParticleFlow/interface/PFCandidate.h +++ b/DataFormats/L1TParticleFlow/interface/PFCandidate.h @@ -48,9 +48,14 @@ namespace l1t { void setZ0(float z0) { setVertex(reco::Particle::Point(0, 0, z0)); } void setDxy(float dxy) { dxy_ = dxy; } + void setCaloEta(float caloeta) { caloEta_ = caloeta; } + void setCaloPhi(float calophi) { caloPhi_ = calophi; } + float z0() const { return vz(); } float dxy() const { return dxy_; } + float caloEta() const { return caloEta_; } + float caloPhi() const { return caloPhi_; } int16_t hwZ0() const { return hwZ0_; } int16_t hwDxy() const { return hwDxy_; } @@ -70,7 +75,7 @@ namespace l1t { PFClusterRef clusterRef_; PFTrackRef trackRef_; MuonRef muonRef_; - float dxy_, puppiWeight_; + float dxy_, puppiWeight_, caloEta_, caloPhi_; int16_t hwZ0_, hwDxy_; uint16_t hwTkQuality_, hwPuppiWeight_, hwEmID_; diff --git a/DataFormats/L1TParticleFlow/interface/layer1_emulator.h b/DataFormats/L1TParticleFlow/interface/layer1_emulator.h index 8262d6b645940..a63e7df8c2a5f 100644 --- a/DataFormats/L1TParticleFlow/interface/layer1_emulator.h +++ b/DataFormats/L1TParticleFlow/interface/layer1_emulator.h @@ -196,6 +196,7 @@ namespace l1ct { const l1t::PFTrack *srcTrack = nullptr; // we use an index to the standalone object needed to retrieve a Ref when putting int sta_idx; + float bdtScore; bool read(std::fstream &from); bool write(std::fstream &to) const; void clear() { @@ -203,6 +204,7 @@ namespace l1ct { srcCluster = nullptr; srcTrack = nullptr; sta_idx = -1; + bdtScore = -999; clearIsoVars(); } diff --git a/DataFormats/L1TParticleFlow/src/PFCandidate.cc b/DataFormats/L1TParticleFlow/src/PFCandidate.cc index 964e55687ec86..36907c32df2ad 100644 --- a/DataFormats/L1TParticleFlow/src/PFCandidate.cc +++ b/DataFormats/L1TParticleFlow/src/PFCandidate.cc @@ -5,6 +5,8 @@ l1t::PFCandidate::PFCandidate( : L1Candidate(p, hwpt, hweta, hwphi, /*hwQuality=*/int(kind)), dxy_(0), puppiWeight_(puppiWeight), + caloEta_(0), + caloPhi_(0), hwZ0_(0), hwDxy_(0), hwTkQuality_(0), diff --git a/DataFormats/L1TParticleFlow/src/classes_def.xml b/DataFormats/L1TParticleFlow/src/classes_def.xml index 5f90ab44bb63d..ff532e682e559 100644 --- a/DataFormats/L1TParticleFlow/src/classes_def.xml +++ b/DataFormats/L1TParticleFlow/src/classes_def.xml @@ -19,7 +19,8 @@ - + + @@ -67,4 +68,3 @@ - diff --git a/L1Trigger/L1THGCalUtilities/test/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc b/L1Trigger/L1THGCalUtilities/test/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc index 34c72aa4df886..6a495b3e5e96e 100644 --- a/L1Trigger/L1THGCalUtilities/test/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc +++ b/L1Trigger/L1THGCalUtilities/test/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc @@ -181,12 +181,12 @@ void HGCalTriggerNtupleHGCMulticlusters::fill(const edm::Event& e, const HGCalTr cl3d_layer_pt_.emplace_back(layer_pt); } - // Retrieve indices of trigger cells inside cluster - cl3d_clusters_id_.emplace_back(cl3d_itr->constituents().size()); - std::transform(cl3d_itr->constituents_begin(), - cl3d_itr->constituents_end(), - cl3d_clusters_id_.back().begin(), - [](const std::pair>& id_cl) { return id_cl.second->detId(); }); + // // Retrieve indices of trigger cells inside cluster + // cl3d_clusters_id_.emplace_back(cl3d_itr->constituents().size()); + // std::transform(cl3d_itr->constituents_begin(), + // cl3d_itr->constituents_end(), + // cl3d_clusters_id_.back().begin(), + // [](const std::pair>& id_cl) { return id_cl.second->detId(); }); } } diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h index b75be595ee9df..93a86c5f4fb8d 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h @@ -4,6 +4,7 @@ #include "DataFormats/L1TParticleFlow/interface/layer1_emulator.h" #include "DataFormats/L1TParticleFlow/interface/egamma.h" #include "DataFormats/L1TParticleFlow/interface/pf.h" +#include "L1Trigger/Phase2L1ParticleFlow/interface/conifer.h" namespace edm { class ParameterSet; @@ -53,6 +54,48 @@ namespace l1ct { bool doPfIso; EGIsoEleObjEmu::IsoType hwIsoTypeTkEle; EGIsoObjEmu::IsoType hwIsoTypeTkEm; + + bool doCompositeTkEle; + struct CompIDParameters { + CompIDParameters(const edm::ParameterSet &); + CompIDParameters(double hoeMin, double hoeMax, double tkptMin, double tkptMax, double srrtotMin, double srrtotMax, double detaMin, double detaMax, double dptMin, double dptMax, double meanzMin, double meanzMax, double dphiMin, double dphiMax, double tkchi2Min, double tkchi2Max, double tkz0Min, double tkz0Max, double tknstubsMin, double tknstubsMax, double BDTcut_wp97p5, double BDTcut_wp95p0) + : hoeMin(hoeMin), hoeMax(hoeMax), + tkptMin(tkptMin),tkptMax(tkptMax), + srrtotMin(srrtotMin),srrtotMax(srrtotMax), + detaMin(detaMin),detaMax(detaMax), + dptMin(dptMin),dptMax(dptMax), + meanzMin(meanzMin),meanzMax(meanzMax), + dphiMin(dphiMin),dphiMax(dphiMax), + tkchi2Min(tkchi2Min),tkchi2Max(tkchi2Max), + tkz0Min(tkz0Min),tkz0Max(tkz0Max), + tknstubsMin(tknstubsMin),tknstubsMax(tknstubsMax), + BDTcut_wp97p5(BDTcut_wp97p5),BDTcut_wp95p0(BDTcut_wp95p0){} + double hoeMin; + double hoeMax; + double tkptMin; + double tkptMax; + double srrtotMin; + double srrtotMax; + double detaMin; + double detaMax; + double dptMin; + double dptMax; + double meanzMin; + double meanzMax; + double dphiMin; + double dphiMax; + double tkchi2Min; + double tkchi2Max; + double tkz0Min; + double tkz0Max; + double tknstubsMin; + double tknstubsMax; + double BDTcut_wp97p5; + double BDTcut_wp95p0; + }; + + CompIDParameters myCompIDparams; + int debug = 0; PFTkEGAlgoEmuConfig(const edm::ParameterSet &iConfig); @@ -81,6 +124,8 @@ namespace l1ct { bool doPfIso = false, EGIsoEleObjEmu::IsoType hwIsoTypeTkEle = EGIsoEleObjEmu::IsoType::TkIso, EGIsoObjEmu::IsoType hwIsoTypeTkEm = EGIsoObjEmu::IsoType::TkIsoPV, + bool doCompositeTkEle = false, + const CompIDParameters &myCompIDparams = {-1.0, 1566.547607421875, 1.9501149654388428, 11102.0048828125, 0.0, 0.01274710614234209, -0.24224889278411865, 0.23079538345336914, 0.010325592942535877, 184.92538452148438, 325.0653991699219, 499.6089782714844, -6.281332015991211, 6.280326843261719, 0.024048099294304848, 1258.37158203125, -14.94140625, 14.94140625, 4.0, 6.0, 0.5406244, 0.9693441}, int debug = 0) : nTRACK(nTrack), @@ -95,6 +140,9 @@ namespace l1ct { emClusterPtMin(emClusterPtMin), dEtaMaxBrem(dEtaMaxBrem), dPhiMaxBrem(dPhiMaxBrem), + //absEtaBoundaries(std::move(absEtaBoundaries)), + //dEtaValues(std::move(dEtaValues)), + //dPhiValues(std::move(dPhiValues)), absEtaBoundaries(absEtaBoundaries), dEtaValues(dEtaValues), dPhiValues(dPhiValues), @@ -108,12 +156,16 @@ namespace l1ct { doPfIso(doPfIso), hwIsoTypeTkEle(hwIsoTypeTkEle), hwIsoTypeTkEm(hwIsoTypeTkEm), + doCompositeTkEle(doCompositeTkEle), + myCompIDparams(myCompIDparams), debug(debug) {} }; class PFTkEGAlgoEmulator { public: - PFTkEGAlgoEmulator(const PFTkEGAlgoEmuConfig &config) : cfg(config), debug_(cfg.debug) {} + + PFTkEGAlgoEmulator(const PFTkEGAlgoEmuConfig &config); + virtual ~PFTkEGAlgoEmulator() {} @@ -133,12 +185,31 @@ namespace l1ct { bool writeEgSta() const { return cfg.writeEgSta; } private: + + void link_emCalo2emCalo(const std::vector &emcalo, std::vector &emCalo2emCalo) const; - void link_emCalo2tk(const PFRegionEmu &r, - const std::vector &emcalo, - const std::vector &track, - std::vector &emCalo2tk) const; + void link_emCalo2tk_elliptic(const PFRegionEmu &r, + const std::vector &emcalo, + const std::vector &track, + std::vector &emCalo2tk) const; + + void link_emCalo2tk_composite(const PFRegionEmu &r, + const std::vector &emcalo, + const std::vector &track, + std::vector &emCalo2tk, + std::vector &emCaloTkBdtScore) const; + + struct CompositeCandidate { + unsigned int cluster_idx; + unsigned int track_idx; + double dpt; // For sorting + }; + + float compute_composite_score(CompositeCandidate &cand, + const std::vector &emcalo, + const std::vector &track, + const PFTkEGAlgoEmuConfig::CompIDParameters ¶ms) const; //FIXME: still needed float deltaPhi(float phi1, float phi2) const; @@ -152,6 +223,7 @@ namespace l1ct { const std::vector &track, const std::vector &emCalo2emCalo, const std::vector &emCalo2tk, + const std::vector &emCaloTkBdtScore, std::vector &egstas, std::vector &egobjs, std::vector &egeleobjs) const; @@ -165,6 +237,7 @@ namespace l1ct { const unsigned int hwQual, const pt_t ptCorr, const int tk_idx, + const float bdtScore, const std::vector &components = {}) const; EGObjEmu &addEGStaToPF(std::vector &egobjs, @@ -182,7 +255,8 @@ namespace l1ct { const EmCaloObjEmu &calo, const TkObjEmu &track, const unsigned int hwQual, - const pt_t ptCorr) const; + const pt_t ptCorr, + const float bdtScore) const; // FIXME: reimplemented from PFAlgoEmulatorBase template @@ -309,6 +383,7 @@ namespace l1ct { z0_t z0) const; PFTkEGAlgoEmuConfig cfg; + std::unique_ptr,ap_fixed<22,3,AP_RND_CONV,AP_SAT>,0>> composite_bdt_; int debug_; }; } // namespace l1ct diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc index f5599d121d0b9..9cade17b5c882 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc @@ -114,6 +114,7 @@ class L1TCorrelatorLayer1Producer : public edm::stream::EDProducer<> { std::unique_ptr fetchEmCalo() const; std::unique_ptr fetchTracks() const; std::unique_ptr fetchPF() const; + std::unique_ptr> fetchDecodedTracks() const; void putPuppi(edm::Event &iEvent) const; void putEgStaObjects(edm::Event &iEvent, @@ -175,6 +176,9 @@ L1TCorrelatorLayer1Producer::L1TCorrelatorLayer1Producer(const edm::ParameterSet #if 0 // LATER produces("TKVtx"); #endif +#if 1 // LATER + produces>("DecodedTK"); +#endif for (const auto &tag : iConfig.getParameter>("emClusters")) { emCands_.push_back(consumes(tag)); @@ -364,6 +368,10 @@ void L1TCorrelatorLayer1Producer::produce(edm::Event &iEvent, const edm::EventSe iEvent.put(fetchEmCalo(), "EmCalo"); iEvent.put(fetchHadCalo(), "Calo"); iEvent.put(fetchTracks(), "TK"); + + #if 1 + iEvent.put(fetchDecodedTracks(), "DecodedTK"); + #endif // Then do the vertexing, and save it out std::vector z0s; @@ -841,6 +849,37 @@ std::unique_ptr L1TCorrelatorLayer1Producer::fetchTr return ret; } +std::unique_ptr> L1TCorrelatorLayer1Producer::fetchDecodedTracks() const { + auto ret = std::make_unique>(); + for (const auto r : event_.decoded.track) { + const auto ® = r.region; + for (const auto &p : r.obj) { + if (p.hwPt == 0 || !reg.isFiducial(p)) + continue; + reco::Particle::PolarLorentzVector p4(p.floatPt(), reg.floatGlbEta(p.hwVtxEta()), reg.floatGlbPhi(p.hwVtxPhi()), 0); + + reco::Particle::Point vtx(0,0,p.floatZ0()); + + ret->emplace_back(l1t::PFTrack(p.intCharge(), + reco::Particle::LorentzVector(p4), + vtx, + p.src->track(), + 0, + reg.floatGlbEta(p.hwEta), + reg.floatGlbPhi(p.hwPhi), + -1, + -1, + p.hwQuality.to_int(), + false, + p.intPt(), + p.intEta(), + p.intPhi())); + } + } + return ret; +} + + std::unique_ptr L1TCorrelatorLayer1Producer::fetchPF() const { auto ret = std::make_unique(); for (unsigned int ir = 0, nr = event_.pfinputs.size(); ir < nr; ++ir) { @@ -861,6 +900,9 @@ std::unique_ptr L1TCorrelatorLayer1Producer::fetchPF ret->back().setHwZ0(p.hwZ0); ret->back().setHwDxy(p.hwDxy); ret->back().setHwTkQuality(p.hwTkQuality); + ret->back().setCaloEta(reg.floatGlbEtaOf(p)); + ret->back().setCaloPhi(reg.floatGlbPhiOf(p)); + setRefs_(ret->back(), p); } for (const auto &p : event_.out[ir].pfneutral) { @@ -871,6 +913,8 @@ std::unique_ptr L1TCorrelatorLayer1Producer::fetchPF p.hwId.isPhoton() ? l1t::PFCandidate::Photon : l1t::PFCandidate::NeutralHadron; ret->emplace_back(type, 0, p4, 1, p.intPt(), p.intEta(), p.intPhi()); ret->back().setHwEmID(p.hwEmID); + ret->back().setCaloEta(reg.floatGlbEtaOf(p)); + ret->back().setCaloPhi(reg.floatGlbPhiOf(p)); setRefs_(ret->back(), p); } } @@ -1042,6 +1086,7 @@ void L1TCorrelatorLayer1Producer::putEgObjects(edm::Event &iEvent, tkele.setHwQual(egele.hwQual); tkele.setPFIsol(egele.floatRelIso(l1ct::EGIsoEleObjEmu::IsoType::PfIso)); tkele.setEgBinaryWord(egele.pack()); + tkele.setCompositeBdtScore(egele.bdtScore); tkeles->push_back(tkele); nele_obj.push_back(tkeles->size() - 1); } diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py index 3d1621b655d56..922fc245e81d2 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py @@ -50,7 +50,32 @@ doTkIso=cms.bool(True), doPfIso=cms.bool(True), hwIsoTypeTkEle=cms.uint32(0), - hwIsoTypeTkEm=cms.uint32(2) + hwIsoTypeTkEm=cms.uint32(2), + doCompositeTkEle=cms.bool(False), + compositeParametersTkEle=cms.PSet( # Parameters used to normalize input features + hoeMin=cms.double(-1.0), + hoeMax=cms.double(1566.547607421875), + tkptMin=cms.double(1.9501149654388428), + tkptMax=cms.double(11102.0048828125), + srrtotMin=cms.double(0.0), + srrtotMax=cms.double(0.01274710614234209), + detaMin=cms.double(-0.24224889278411865), + detaMax=cms.double(0.23079538345336914), + dptMin=cms.double(0.010325592942535877), + dptMax=cms.double(184.92538452148438), + meanzMin=cms.double(325.0653991699219), + meanzMax=cms.double(499.6089782714844), + dphiMin=cms.double(-6.281332015991211), + dphiMax=cms.double(6.280326843261719), + tkchi2Min=cms.double(0.024048099294304848), + tkchi2Max=cms.double(1258.37158203125), + tkz0Min=cms.double(-14.94140625), + tkz0Max=cms.double(14.94140625), + tknstubsMin=cms.double(4.0), + tknstubsMax=cms.double(6.0), + BDTcut_wp97p5=cms.double(0.5406244), + BDTcut_wp95p0=cms.double(0.9693441), + ), ) tkEgSorterParameters = cms.PSet( diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py index e0b2e6542ae5d..2464213e50755 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py @@ -261,7 +261,8 @@ doBremRecovery=True, doEndcapHwQual=True, writeBeforeBremRecovery=False, - writeEGSta=True), + writeEGSta=True, + doCompositeTkEle=False), tkEgSorterParameters=tkEgSorterParameters.clone( nObjToSort = 5 ), diff --git a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp index 2fc881879455d..9773754ced2e7 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp @@ -8,6 +8,8 @@ #include #include +#include "L1Trigger/Phase2L1ParticleFlow/src/dbgPrintf.h" +#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" using namespace l1ct; #ifdef CMSSW_GIT_HASH @@ -39,6 +41,8 @@ l1ct::PFTkEGAlgoEmuConfig::PFTkEGAlgoEmuConfig(const edm::ParameterSet &pset) doPfIso(pset.getParameter("doPfIso")), hwIsoTypeTkEle(static_cast(pset.getParameter("hwIsoTypeTkEle"))), hwIsoTypeTkEm(static_cast(pset.getParameter("hwIsoTypeTkEm"))), + doCompositeTkEle(pset.getParameter("doCompositeTkEle")), + myCompIDparams(pset.getParameter("compositeParametersTkEle")), debug(pset.getUntrackedParameter("debug", 0)) {} l1ct::PFTkEGAlgoEmuConfig::IsoParameters::IsoParameters(const edm::ParameterSet &pset) @@ -47,8 +51,42 @@ l1ct::PFTkEGAlgoEmuConfig::IsoParameters::IsoParameters(const edm::ParameterSet pset.getParameter("dRMin"), pset.getParameter("dRMax")) {} +l1ct::PFTkEGAlgoEmuConfig::CompIDParameters::CompIDParameters(const edm::ParameterSet &pset) + : CompIDParameters(pset.getParameter("hoeMin"), + pset.getParameter("hoeMax"), + pset.getParameter("tkptMin"), + pset.getParameter("tkptMax"), + pset.getParameter("srrtotMin"), + pset.getParameter("srrtotMax"), + pset.getParameter("detaMin"), + pset.getParameter("detaMax"), + pset.getParameter("dptMin"), + pset.getParameter("dptMax"), + pset.getParameter("meanzMin"), + pset.getParameter("meanzMax"), + pset.getParameter("dphiMin"), + pset.getParameter("dphiMax"), + pset.getParameter("tkchi2Min"), + pset.getParameter("tkchi2Max"), + pset.getParameter("tkz0Min"), + pset.getParameter("tkz0Max"), + pset.getParameter("tknstubsMin"), + pset.getParameter("tknstubsMax"), + pset.getParameter("BDTcut_wp97p5"), + pset.getParameter("BDTcut_wp95p0")) {} + #endif +PFTkEGAlgoEmulator::PFTkEGAlgoEmulator(const PFTkEGAlgoEmuConfig &config) : cfg(config), +composite_bdt_(nullptr), +debug_(cfg.debug) { + if(cfg.doCompositeTkEle) { + //FIXME: make the name of the file configurable + auto resolvedFileName = edm::FileInPath("L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/compositeID.json").fullPath(); + composite_bdt_ = std::make_unique,ap_fixed<22,3,AP_RND_CONV,AP_SAT>,0>> (resolvedFileName); + } +} + void PFTkEGAlgoEmulator::toFirmware(const PFInputRegion &in, PFRegion ®ion, EmCaloObj emcalo[/*nCALO*/], @@ -109,10 +147,11 @@ void PFTkEGAlgoEmulator::link_emCalo2emCalo(const std::vector &emc } } -void PFTkEGAlgoEmulator::link_emCalo2tk(const PFRegionEmu &r, - const std::vector &emcalo, - const std::vector &track, - std::vector &emCalo2tk) const { + +void PFTkEGAlgoEmulator::link_emCalo2tk_elliptic(const PFRegionEmu &r, + const std::vector &emcalo, + const std::vector &track, + std::vector &emCalo2tk) const { unsigned int nTrackMax = std::min(track.size(), cfg.nTRACK_EGIN); for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) { auto &calo = emcalo[ic]; @@ -146,6 +185,117 @@ void PFTkEGAlgoEmulator::link_emCalo2tk(const PFRegionEmu &r, } } + +void PFTkEGAlgoEmulator::link_emCalo2tk_composite(const PFRegionEmu &r, + const std::vector &emcalo, + const std::vector &track, + std::vector &emCalo2tk, + std::vector &emCaloTkBdtScore) const { + //FIXME: should be configurable + const int nCAND_PER_CLUSTER = 4; + unsigned int nTrackMax = std::min(track.size(), cfg.nTRACK_EGIN); + for (int ic = 0, nc = emcalo.size(); ic < nc; ++ic) { + auto &calo = emcalo[ic]; + + std::vector candidates; + + for (unsigned int itk = 0; itk < nTrackMax; ++itk) { + const auto &tk = track[itk]; + if (tk.floatPt() < cfg.trkQualityPtMin) + continue; + + float d_phi = deltaPhi(tk.floatPhi(), calo.floatPhi()); + float d_eta = tk.floatEta() - calo.floatEta(); // We only use it squared + float dR = sqrt((d_phi * d_phi ) + (d_eta * d_eta )); + + if (dR<0.2){ + // Only store indices, dR and dpT for now. The other quantities are computed only for the best nCandPerCluster. + CompositeCandidate cand; + cand.cluster_idx = ic; + cand.track_idx = itk; + cand.dpt = fabs(tk.floatPt() - calo.floatPt()); + candidates.push_back(cand); + } + } + // FIXME: find best sort criteria, for now we use dpt + std::sort(candidates.begin(), candidates.end(), + [](const CompositeCandidate & a, const CompositeCandidate & b) -> bool + { return a.dpt < b.dpt; }); + unsigned int nCandPerCluster = std::min(candidates.size(), nCAND_PER_CLUSTER); + std::cout << "# composit candidates: " << nCandPerCluster << std::endl; + if(nCandPerCluster == 0) continue; + + float bdtWP_MVA = cfg.myCompIDparams.BDTcut_wp97p5; + float bdtWP_XGB = 1. / (1. + std::sqrt((1. - bdtWP_MVA) / (1. + bdtWP_MVA))); // Convert WP value from ROOT.TMVA to XGboost + float maxScore = -999; + int ibest = -1; + for(unsigned int icand = 0; icand < nCandPerCluster; icand++) { + auto &cand = candidates[icand]; + std::vector emcalo_sel = emcalo; + float score = compute_composite_score(cand, emcalo_sel, track, cfg.myCompIDparams); + if(score > maxScore) { + // if((score > bdtWP_XGB) && (score > maxScore)) { + maxScore = score; + ibest = icand; + } + } + if(ibest != -1) { + emCalo2tk[ic] = candidates[ibest].track_idx; + emCaloTkBdtScore[ic] = maxScore; + } + } +} + + +float PFTkEGAlgoEmulator::compute_composite_score(CompositeCandidate &cand, + const std::vector &emcalo, + const std::vector &track, + const PFTkEGAlgoEmuConfig::CompIDParameters ¶ms) const { + // Get the cluster/track objects that form the composite candidate + const auto &calo = emcalo[cand.cluster_idx]; + const auto &tk = track[cand.track_idx]; + + // FIXME: using these two floats crashes code... + float srrtot_f = dynamic_cast(calo.src->constituentsAndFractions()[0].first.get())->sigmaRRTot(); + float meanz_f = abs(dynamic_cast(calo.src->constituentsAndFractions()[0].first.get())->zBarycenter()); + + // Call and normalize input feature values, then cast to ap_fixed. + // Note that for some features (e.g. track pT) we call the floating point representation, but that's already quantized! + // Several other features, such as chi2 or most cluster features, are not quantized before casting them to ap_fixed. + ap_fixed<22,3,AP_RND_CONV,AP_SAT> hoe = (calo.src->hOverE()-params.hoeMin)/(params.hoeMax-params.hoeMin); + ap_fixed<22,3,AP_RND_CONV,AP_SAT> tkpt = (tk.floatPt()-params.tkptMin)/(params.tkptMax-params.tkptMin); + ap_fixed<22,3,AP_RND_CONV,AP_SAT> srrtot = (srrtot_f-params.srrtotMin)/(params.srrtotMax-params.srrtotMin); + ap_fixed<22,3,AP_RND_CONV,AP_SAT> deta = (tk.floatEta() - calo.floatEta()-params.detaMin)/(params.detaMax-params.detaMin); + // FIXME: do we really need dpt to be a ratio? + ap_fixed<22,3,AP_RND_CONV,AP_SAT> dpt = ((tk.floatPt()/calo.floatPt())-params.dptMin)/(params.dptMax-params.dptMin); + ap_fixed<22,3,AP_RND_CONV,AP_SAT> meanz = (meanz_f-params.meanzMin)/(params.meanzMax-params.meanzMin); + ap_fixed<22,3,AP_RND_CONV,AP_SAT> dphi = (deltaPhi(tk.floatPhi(), calo.floatPhi()) -params.dphiMin)/(params.dphiMax-params.dphiMin); + ap_fixed<22,3,AP_RND_CONV,AP_SAT> chi2 = (tk.src->chi2()-params.tkchi2Min)/(params.tkchi2Max-params.tkchi2Min); + ap_fixed<22,3,AP_RND_CONV,AP_SAT> tkz0 = (tk.floatZ0()-params.tkz0Min)/(params.tkz0Max-params.tkz0Min); + ap_fixed<22,3,AP_RND_CONV,AP_SAT> nstubs = (tk.src->nStubs()-params.tknstubsMin)/(params.tknstubsMax-params.tknstubsMin); + // std::cout<<"hoe\t"<hOverE()<<"\t"<<(calo.src->hOverE()-params.hoeMin)/(params.hoeMax-params.hoeMin)<chi2()<<"\t"<<(tk.src->chi2()-params.tkchi2Min)/(params.tkchi2Max-params.tkchi2Min)<nStubs()<<"\t"<<(tk.src->nStubs()-params.tknstubsMin)/(params.tknstubsMax-params.tknstubsMin)<> inputs = { hoe, tkpt, srrtot, deta, dpt, meanz, dphi, chi2, tkz0, nstubs } ; + auto bdt_score = composite_bdt_->decision_function(inputs); + + float bdt_score_CON = bdt_score[0]; + float bdt_score_XGB = 1/(1+exp(-bdt_score_CON)); // Map Conifer score to XGboost score. (same as scipy.expit) + + // std::cout<<"BDT score of composite candidate = "< &emcalo, std::vector &emcalo_sel) const { @@ -183,12 +333,18 @@ void PFTkEGAlgoEmulator::run(const PFInputRegion &in, OutputRegion &out) const { link_emCalo2emCalo(emcalo_sel, emCalo2emCalo); std::vector emCalo2tk(emcalo_sel.size(), -1); - link_emCalo2tk(in.region, emcalo_sel, in.track, emCalo2tk); + std::vector emCaloTkBdtScore(emcalo_sel.size(), -999); + if(cfg.doCompositeTkEle) { + link_emCalo2tk_composite(in.region, emcalo_sel, in.track, emCalo2tk, emCaloTkBdtScore); + } else { + link_emCalo2tk_elliptic(in.region, emcalo_sel, in.track, emCalo2tk); + } + out.egsta.clear(); std::vector egobjs; std::vector egeleobjs; - eg_algo(in.region, emcalo_sel, in.track, emCalo2emCalo, emCalo2tk, out.egsta, egobjs, egeleobjs); + eg_algo(in.region, emcalo_sel, in.track, emCalo2emCalo, emCalo2tk, emCaloTkBdtScore, out.egsta, egobjs, egeleobjs); unsigned int nEGOut = std::min(cfg.nEM_EGOUT, egobjs.size()); unsigned int nEGEleOut = std::min(cfg.nEM_EGOUT, egeleobjs.size()); @@ -205,6 +361,7 @@ void PFTkEGAlgoEmulator::eg_algo(const PFRegionEmu ®ion, const std::vector &track, const std::vector &emCalo2emCalo, const std::vector &emCalo2tk, + const std::vector &emCaloTkBdtScore, std::vector &egstas, std::vector &egobjs, std::vector &egeleobjs) const { @@ -220,6 +377,7 @@ void PFTkEGAlgoEmulator::eg_algo(const PFRegionEmu ®ion, << " phi " << calo.hwPhi << std::endl; int itk = emCalo2tk[ic]; + float bdt = emCaloTkBdtScore[ic]; // check if brem recovery is on if (!cfg.doBremRecovery || cfg.writeBeforeBremRecovery) { @@ -232,7 +390,7 @@ void PFTkEGAlgoEmulator::eg_algo(const PFRegionEmu ®ion, egQual = calo.hwEmID | 0x8; } - addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, egQual, calo.hwPt, itk); + addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, egQual, calo.hwPt, itk, bdt); } if (!cfg.doBremRecovery) @@ -255,7 +413,7 @@ void PFTkEGAlgoEmulator::eg_algo(const PFRegionEmu ®ion, // 2. create EG objects with brem recovery // NOTE: duplicating the object is suboptimal but this is done for keeping things as in TDR code... - addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, calo.hwEmID, ptBremReco, itk, components); + addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, calo.hwEmID, ptBremReco, itk, bdt, components); } } @@ -309,7 +467,8 @@ EGIsoEleObjEmu &PFTkEGAlgoEmulator::addEGIsoEleToPF(std::vector const EmCaloObjEmu &calo, const TkObjEmu &track, const unsigned int hwQual, - const pt_t ptCorr) const { + const pt_t ptCorr, + const float bdtScore) const { EGIsoEleObjEmu egiso; egiso.clear(); egiso.hwPt = ptCorr; @@ -328,6 +487,7 @@ EGIsoEleObjEmu &PFTkEGAlgoEmulator::addEGIsoEleToPF(std::vector egiso.hwCharge = track.hwCharge; egiso.srcCluster = calo.src; egiso.srcTrack = track.src; + egiso.bdtScore = bdtScore; egobjs.push_back(egiso); if (debug_ > 2) @@ -346,6 +506,7 @@ void PFTkEGAlgoEmulator::addEgObjsToPF(std::vector &egstas, const unsigned int hwQual, const pt_t ptCorr, const int tk_idx, + const float bdtScore, const std::vector &components) const { int sta_idx = -1; if (writeEgSta()) { @@ -355,7 +516,7 @@ void PFTkEGAlgoEmulator::addEgObjsToPF(std::vector &egstas, EGIsoObjEmu &egobj = addEGIsoToPF(egobjs, emcalo[calo_idx], hwQual, ptCorr); egobj.sta_idx = sta_idx; if (tk_idx != -1) { - EGIsoEleObjEmu &eleobj = addEGIsoEleToPF(egeleobjs, emcalo[calo_idx], track[tk_idx], hwQual, ptCorr); + EGIsoEleObjEmu &eleobj = addEGIsoEleToPF(egeleobjs, emcalo[calo_idx], track[tk_idx], hwQual, ptCorr, bdtScore); eleobj.sta_idx = sta_idx; } } diff --git a/L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/compositeID.json b/L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/compositeID.json new file mode 100644 index 0000000000000..c7cfb823abdd1 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/compositeID.json @@ -0,0 +1 @@ +{"max_depth": 3, "n_trees": 10, "n_classes": 2, "n_features": 10, "trees": [[{"feature": [0, 2, 0, 1, 1, 1, 0, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000703328056, 0.332558006, 0.000770003942, 0.000253851322, 0.000915614248, 0.000717727817, 0.000839613494, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.0352521278, 0.590574145, -0.54147464, 0.185335979, -0.503847778, 0.257735074, -0.503176689, -0.595144987], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [0, 2, 0, 1, 1, 1, 5, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000697391108, 0.340182066, 0.000747981714, 0.000371906266, 0.00168169127, 0.000789736281, 0.0738627762, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.0733107328, 0.458879024, -0.396272093, 0.309196174, -0.360979736, 0.275205106, -0.349665552, -0.460340261], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [0, 1, 1, 4, 2, 5, 5, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000694462447, 0.000501940958, 0.00114817475, 0.000432471978, 0.370490968, 0.0738627762, 0.0940087885, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.380878597, -0.283362597, 0.397695929, -0.294095337, -0.280492604, -0.399287552, 0.310794801, -0.362855673], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [0, 2, 0, 3, 1, 2, 0, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000709932996, 0.321239412, 0.000803797971, 0.602574646, 0.00082699745, 0.336098373, 0.000910258619, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.353067011, -0.383320004, -0.345221847, 0.114074722, 0.156569079, -0.339963675, -0.309075147, -0.367914289], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [0, 1, 1, 4, 2, 7, 2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000687893422, 0.000643459614, 0.000599945546, 0.000337057747, 0.358681321, 0.0250910092, 0.343107432, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.349036545, -0.16346851, 0.341645479, -0.123347722, -0.342872322, -0.239095137, 0.228956148, -0.307008922], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [0, 1, 0, 4, 2, 1, 0, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000722329598, 0.000350446324, 0.000839613494, 0.000351717492, 0.365002632, 0.00043548242, 0.000917885685, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.253911734, -0.244474351, 0.316286922, -0.238444999, -0.279888034, 0.0941597223, -0.259838462, -0.331297874], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [0, 3, 5, 8, 3, 1, 1, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000677733915, 0.438798308, 0.0826261193, 0.850980401, 0.592054725, 0.00032415273, 0.00283236895, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, -0.413880289, 0.257492244, 0.303176314, -0.260525465, -0.260994583, 0.136501729, -0.314218611, 0.153914765], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [1, 0, 2, 2, 0, 5, 2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.00122120825, 0.000803797971, 0.33025986, 0.301885217, 0.000917885685, 0.133645415, 0.387662947, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.145611659, -0.207396939, -0.205260798, -0.312455952, 0.311472386, -0.278300524, -0.0131324222, -0.31602335], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [1, 5, 2, 2, 5, 0, 4, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.00148988748, 0.0920395404, 0.357989103, 0.310617983, 0.116608188, 0.00109421043, 0.00491232192, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.13668026, -0.192945674, -0.225791544, -0.311047494, 0.302013129, -0.30914408, -0.0594445392, -0.352083832], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [1, 7, 2, 4, 0, 5, 2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000719727308, 0.0146044344, 0.315754294, 0.000207430043, 0.000914971635, 0.114892721, 0.404698849, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.209280849, -0.279667407, 0.159891739, -0.288006365, 0.288983703, -0.244072288, 0.0150486836, -0.272770345], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}]], "init_predict": [0, 0], "norm": 1} \ No newline at end of file diff --git a/L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/conifer.h b/L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/conifer.h new file mode 100644 index 0000000000000..d9529aa9dcb90 --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/conifer.h @@ -0,0 +1,158 @@ +#ifndef CONIFER_CPP_H__ +#define CONIFER_CPP_H__ +#include "nlohmann/json.hpp" +#include +#include + +namespace conifer{ + +/* --- +* Balanced tree reduce implementation. +* Reduces an array of inputs to a single value using the template binary operator 'Op', +* for example summing all elements with Op_add, or finding the maximum with Op_max +* Use only when the input array is fully unrolled. Or, slice out a fully unrolled section +* before applying and accumulate the result over the rolled dimension. +* Required for emulation to guarantee equality of ordering. +* --- */ +constexpr int floorlog2(int x) { return (x < 2) ? 0 : 1 + floorlog2(x / 2); } + +template +constexpr int pow(int x) { + return x == 0 ? 1 : B * pow(x - 1); +} + +constexpr int pow2(int x) { return pow<2>(x); } + +template +T reduce(std::vector x, Op op) { + int N = x.size(); + int leftN = pow2(floorlog2(N - 1)) > 0 ? pow2(floorlog2(N - 1)) : 0; + //static constexpr int rightN = N - leftN > 0 ? N - leftN : 0; + if (N == 1) { + return x.at(0); + } else if (N == 2) { + return op(x.at(0), x.at(1)); + } else { + std::vector left(x.begin(), x.begin() + leftN); + std::vector right(x.begin() + leftN, x.end()); + return op(reduce(left, op), reduce(right, op)); + } +} + +template +class OpAdd { +public: + T operator()(T a, T b) { return a + b; } +}; + +template +class DecisionTree{ + +private: + std::vector feature; + std::vector children_left; + std::vector children_right; + std::vector threshold_; + std::vector value_; + std::vector threshold; + std::vector value; + +public: + + U decision_function(std::vector x) const{ + /* Do the prediction */ + int i = 0; + while(feature[i] != -2){ // continue until reaching leaf + bool comparison = x[feature[i]] <= threshold_[i]; + i = comparison ? children_left[i] : children_right[i]; + } + return value_[i]; + } + + void init_(){ + /* Since T, U types may not be readable from the JSON, read them to double and the cast them here */ + std::transform(threshold.begin(), threshold.end(), std::back_inserter(threshold_), + [](double t) -> T { return (T) t; }); + std::transform(value.begin(), value.end(), std::back_inserter(value_), + [](double v) -> U { return (U) v; }); + } + + // Define how to read this class to/from JSON + NLOHMANN_DEFINE_TYPE_INTRUSIVE(DecisionTree, feature, children_left, children_right, threshold, value); + +}; // class DecisionTree + +template +class BDT{ + +private: + + int n_classes; + int n_trees; + unsigned int n_features; + std::vector init_predict; + std::vector init_predict_; + // vector of decision trees: outer dimension tree, inner dimension class + std::vector>> trees; + OpAdd add; + +public: + + // Define how to read this class to/from JSON + NLOHMANN_DEFINE_TYPE_INTRUSIVE(BDT, n_classes, n_trees, n_features, init_predict, trees); + + BDT(std::string filename){ + /* Construct the BDT from conifer cpp backend JSON file */ + std::ifstream ifs(filename); + nlohmann::json j = nlohmann::json::parse(ifs); + from_json(j, *this); + /* Do some transformation to initialise things into the proper emulation T, U types */ + if(n_classes == 2) n_classes = 1; + std::transform(init_predict.begin(), init_predict.end(), std::back_inserter(init_predict_), + [](double ip) -> U { return (U) ip; }); + for(int i = 0; i < n_trees; i++){ + for(int j = 0; j < n_classes; j++){ + trees.at(i).at(j).init_(); + } + } + } + + std::vector decision_function(std::vector x) const{ + /* Do the prediction */ + // std::cout< values; + std::vector> values_trees; + values_trees.resize(n_classes); + values.resize(n_classes, U(0)); + for(int i = 0; i < n_classes; i++){ + std::transform(trees.begin(), trees.end(), std::back_inserter(values_trees.at(i)), + [&i, &x](auto tree_v){ return tree_v.at(i).decision_function(x); }); + if(useAddTree){ + values.at(i) = init_predict_.at(i); + values.at(i) += reduce>(values_trees.at(i), add); + }else{ + values.at(i) = std::accumulate(values_trees.at(i).begin(), values_trees.at(i).end(), U(init_predict_.at(i))); + } + } + + return values; + } + + std::vector _decision_function_double(std::vector x) const{ + /* Do the prediction with data in/out as double, cast to T, U before prediction */ + std::vector xt; + std::transform(x.begin(), x.end(), std::back_inserter(xt), + [](double xi) -> T { return (T) xi; }); + std::vector y = decision_function(xt); + std::vector yd; + std::transform(y.begin(), y.end(), std::back_inserter(yd), + [](U yi) -> double { return (double) yi; }); + return yd; + } + +}; // class BDT + +} // namespace conifer + +#endif \ No newline at end of file diff --git a/L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_dumpFiles_cfg.py b/L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_dumpFiles_cfg.py index c8eddfb3b029c..edc30246e11c4 100644 --- a/L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_dumpFiles_cfg.py +++ b/L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_dumpFiles_cfg.py @@ -7,7 +7,7 @@ process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") process.load("FWCore.MessageLogger.MessageLogger_cfi") process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True), allowUnscheduled = cms.untracked.bool(False) ) -process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(1000)) +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(100)) process.MessageLogger.cerr.FwkReport.reportEvery = 1 process.source = cms.Source("PoolSource", From cc0a805434e97828d8057b078b7a9a11feb0f7ab Mon Sep 17 00:00:00 2001 From: Peter Eduard Meiring Date: Fri, 16 Sep 2022 16:32:08 +0200 Subject: [PATCH 02/13] Remove pt10 cut from composite candidates --- L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py | 3 ++- .../Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp | 8 ++++---- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py index 2464213e50755..686995437e5b3 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py @@ -262,7 +262,8 @@ doEndcapHwQual=True, writeBeforeBremRecovery=False, writeEGSta=True, - doCompositeTkEle=False), + doCompositeTkEle=True + trkQualityPtMin=0.), # This should be 10 GeV when doCompositeTkEle=False tkEgSorterParameters=tkEgSorterParameters.clone( nObjToSort = 5 ), diff --git a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp index 9773754ced2e7..ff7f54cfc7135 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp @@ -82,7 +82,7 @@ composite_bdt_(nullptr), debug_(cfg.debug) { if(cfg.doCompositeTkEle) { //FIXME: make the name of the file configurable - auto resolvedFileName = edm::FileInPath("L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/compositeID.json").fullPath(); + auto resolvedFileName = edm::FileInPath("L1Trigger/Phase2L1ParticleFlow/data/compositeID.json").fullPath(); composite_bdt_ = std::make_unique,ap_fixed<22,3,AP_RND_CONV,AP_SAT>,0>> (resolvedFileName); } } @@ -200,9 +200,9 @@ void PFTkEGAlgoEmulator::link_emCalo2tk_composite(const PFRegionEmu &r, std::vector candidates; for (unsigned int itk = 0; itk < nTrackMax; ++itk) { - const auto &tk = track[itk]; - if (tk.floatPt() < cfg.trkQualityPtMin) - continue; + const auto &tk = track[itk]; + if (tk.floatPt() <= cfg.trkQualityPtMin) + continue; float d_phi = deltaPhi(tk.floatPhi(), calo.floatPhi()); float d_eta = tk.floatEta() - calo.floatEta(); // We only use it squared From 959010a88bbea94ac9e94432a1c4cd96cb060764 Mon Sep 17 00:00:00 2001 From: Gianluca Date: Mon, 10 Oct 2022 23:54:24 +0200 Subject: [PATCH 03/13] Adapt inout data-formats for Composite BDT evaluation. --- .../L1TParticleFlow/interface/PFCluster.h | 17 +- .../L1TParticleFlow/interface/datatypes.h | 17 ++ .../interface/layer1_emulator.h | 8 +- .../L1TParticleFlow/interface/layer1_objs.h | 55 +++++- .../L1TParticleFlow/src/classes_def.xml | 3 +- .../Phase2L1ParticleFlow/interface/conifer.h | 14 +- .../interface/egamma/pftkegalgo_ref.h | 17 +- .../plugins/L1TCorrelatorLayer1Producer.cc | 16 +- .../PFClusterProducerFromHGC3DClusters.cc | 4 + .../python/l1TkEgAlgoEmulator_cfi.py | 1 + .../python/l1ctLayer1_cff.py | 2 +- .../src/egamma/pfeginput_ref.cpp | 4 + .../src/egamma/pftkegalgo_ref.cpp | 67 ++++---- .../src/l1-converters/hgcalinputt_ref.cpp | 2 + .../src/newfirmware/egamma/compositeID.json | 1 - .../src/newfirmware/egamma/conifer.h | 158 ------------------ 16 files changed, 162 insertions(+), 224 deletions(-) delete mode 100644 L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/compositeID.json delete mode 100644 L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/conifer.h diff --git a/DataFormats/L1TParticleFlow/interface/PFCluster.h b/DataFormats/L1TParticleFlow/interface/PFCluster.h index 1851746b51d19..2272ccbe353db 100644 --- a/DataFormats/L1TParticleFlow/interface/PFCluster.h +++ b/DataFormats/L1TParticleFlow/interface/PFCluster.h @@ -22,10 +22,14 @@ namespace l1t { float ptError = 0, int hwpt = 0, int hweta = 0, - int hwphi = 0) + int hwphi = 0, + float absZBarycenter = 0., + float sigmaRR = 0.) : L1Candidate(PolarLorentzVector(pt, eta, phi, 0), hwpt, hweta, hwphi, /*hwQuality=*/isEM ? 1 : 0), hOverE_(hOverE), - ptError_(ptError) { + ptError_(ptError), + absZBarycenter_(absZBarycenter), + sigmaRR_(sigmaRR) { setPdgId(isEM ? 22 : 130); // photon : non-photon(K0) } PFCluster( @@ -37,6 +41,12 @@ namespace l1t { float hOverE() const { return hOverE_; } void setHOverE(float hOverE) { hOverE_ = hOverE; } + void setSigmaRR(float sigmaRR) { sigmaRR_ = sigmaRR; } + float absZBarycenter() const { return absZBarycenter_; } + + void setAbsZBarycenter(float absZBarycenter) { absZBarycenter_ = absZBarycenter; } + float sigmaRR() const { return sigmaRR_; } + float emEt() const { if (hOverE_ == -1) return 0; @@ -68,6 +78,9 @@ namespace l1t { private: float hOverE_, ptError_, egVsPionMVAOut_, egVsPUMVAOut_; + // HGC dedicated quantities (0ed by default) + float absZBarycenter_, sigmaRR_; + ConstituentsAndFractions constituents_; }; diff --git a/DataFormats/L1TParticleFlow/interface/datatypes.h b/DataFormats/L1TParticleFlow/interface/datatypes.h index f56e7b36c2495..719ab378d3205 100644 --- a/DataFormats/L1TParticleFlow/interface/datatypes.h +++ b/DataFormats/L1TParticleFlow/interface/datatypes.h @@ -39,6 +39,13 @@ namespace l1ct { typedef ap_uint<10> em2calo_dr_t; typedef ap_uint<13> tk2calo_dq_t; typedef ap_uint<4> egquality_t; + typedef ap_uint<3> stub_t; + // FIXME: random choice for the various bitwidth below!!!! + typedef ap_ufixed<12, 10, AP_TRN, AP_SAT> chi2_t; + typedef ap_ufixed<10, 1, AP_TRN, AP_SAT> srrtot_t; + typedef ap_uint<8> meanz_t; // mean - SCALE_MEANZ = 320 + typedef ap_ufixed<10, 5, AP_TRN, AP_SAT> hoe_t; + // FIXME: adjust range 10-11bits -> 1/4 - 1/2TeV is probably more than enough for all reasonable use cases typedef ap_ufixed<11, 9, AP_TRN, AP_SAT> iso_t; @@ -149,6 +156,8 @@ namespace l1ct { constexpr float Z0_LSB = 0.05; constexpr float DXY_LSB = 0.05; constexpr float PUPPIW_LSB = 1.0 / 256; + constexpr float MEANZ_SCALE = 320.; + inline float floatPt(pt_t pt) { return pt.to_float(); } inline float floatPt(dpt_t pt) { return pt.to_float(); } inline float floatPt(pt2_t pt2) { return pt2.to_float(); } @@ -164,6 +173,10 @@ namespace l1ct { inline float floatDxy(dxy_t dxy) { return dxy.to_float() * DXY_LSB; } inline float floatPuppiW(puppiWgt_t puppiw) { return puppiw.to_float() * PUPPIW_LSB; } inline float floatIso(iso_t iso) { return iso.to_float(); } + inline float floatChi2(chi2_t chi2) { return chi2.to_float(); } + inline float floatSrrTot(srrtot_t srrtot) { return srrtot.to_float(); }; + inline float floatMeanZ(meanz_t meanz) { return meanz + MEANZ_SCALE; }; + inline float floatHoe(hoe_t hoe) { return hoe.to_float(); }; inline pt_t makePt(int pt) { return ap_ufixed<16, 14>(pt) >> 2; } inline dpt_t makeDPt(int dpt) { return ap_fixed<18, 16>(dpt) >> 2; } @@ -194,6 +207,10 @@ namespace l1ct { inline iso_t makeIso(float iso) { return iso_t(0.25 * round(iso * 4)); } inline int makeDR2FromFloatDR(float dr) { return ceil(dr * dr / ETAPHI_LSB / ETAPHI_LSB); } + inline chi2_t makeChi2(float chi2) { return chi2_t(chi2); } + inline srrtot_t makeSrrTot(float var) { return srrtot_t(var); }; + inline meanz_t makeMeanZ(float var) { return round(var - MEANZ_SCALE); }; + inline hoe_t makeHoe(float var) { return hoe_t(var); }; inline float maxAbsEta() { return ((1 << (eta_t::width - 1)) - 1) * ETAPHI_LSB; } inline float maxAbsPhi() { return ((1 << (phi_t::width - 1)) - 1) * ETAPHI_LSB; } diff --git a/DataFormats/L1TParticleFlow/interface/layer1_emulator.h b/DataFormats/L1TParticleFlow/interface/layer1_emulator.h index a63e7df8c2a5f..152456059b939 100644 --- a/DataFormats/L1TParticleFlow/interface/layer1_emulator.h +++ b/DataFormats/L1TParticleFlow/interface/layer1_emulator.h @@ -39,7 +39,7 @@ namespace l1ct { }; struct TkObjEmu : public TkObj { - uint16_t hwChi2, hwStubs; + // uint16_t hwChi2, hwStubs; float simPt, simCaloEta, simCaloPhi, simVtxEta, simVtxPhi, simZ0, simD0; const l1t::PFTrack *src = nullptr; bool read(std::fstream &from); @@ -47,8 +47,8 @@ namespace l1ct { void clear() { TkObj::clear(); src = nullptr; - hwChi2 = 0; - hwStubs = 0; + // hwChi2 = 0; + // hwStubs = 0; simPt = 0; simCaloEta = 0; simCaloPhi = 0; @@ -336,7 +336,7 @@ namespace l1ct { }; struct Event { - enum { VERSION = 11 }; + enum { VERSION = 12 }; uint32_t run, lumi; uint64_t event; RawInputs raw; diff --git a/DataFormats/L1TParticleFlow/interface/layer1_objs.h b/DataFormats/L1TParticleFlow/interface/layer1_objs.h index 83765fc4b2b10..a5aab55dcbf8f 100644 --- a/DataFormats/L1TParticleFlow/interface/layer1_objs.h +++ b/DataFormats/L1TParticleFlow/interface/layer1_objs.h @@ -12,10 +12,13 @@ namespace l1ct { phi_t hwPhi; // relative to the region center, at calo pt_t hwEmPt; emid_t hwEmID; + srrtot_t hwSrrTot; + meanz_t hwMeanZ; + hoe_t hwHoe; inline bool operator==(const HadCaloObj &other) const { return hwPt == other.hwPt && hwEta == other.hwEta && hwPhi == other.hwPhi && hwEmPt == other.hwEmPt && - hwEmID == other.hwEmID; + hwEmID == other.hwEmID && hwSrrTot == other.hwSrrTot && hwMeanZ == other.hwMeanZ && hwHoe == other.hwHoe; } inline bool operator>(const HadCaloObj &other) const { return hwPt > other.hwPt; } @@ -27,6 +30,9 @@ namespace l1ct { hwPhi = 0; hwEmPt = 0; hwEmID = 0; + hwSrrTot = 0; + hwMeanZ = 0; + hwHoe = 0; } int intPt() const { return Scales::intPt(hwPt); } @@ -37,10 +43,14 @@ namespace l1ct { float floatEmPt() const { return Scales::floatPt(hwEmPt); } float floatEta() const { return Scales::floatEta(hwEta); } float floatPhi() const { return Scales::floatPhi(hwPhi); } + float floatSrrTot() const { return Scales::floatSrrTot(hwSrrTot); }; + float floatMeanZ() const { return Scales::floatMeanZ(hwMeanZ); }; + float floatHoe() const { return Scales::floatHoe(hwHoe); }; bool hwIsEM() const { return hwEmID != 0; } - static const int BITWIDTH = pt_t::width + eta_t::width + phi_t::width + pt_t::width + emid_t::width; + static const int BITWIDTH = pt_t::width + eta_t::width + phi_t::width + pt_t::width + emid_t::width + + srrtot_t::width + meanz_t::width + hoe_t::width; inline ap_uint pack() const { ap_uint ret; unsigned int start = 0; @@ -49,6 +59,9 @@ namespace l1ct { pack_into_bits(ret, start, hwPhi); pack_into_bits(ret, start, hwEmPt); pack_into_bits(ret, start, hwEmID); + pack_into_bits(ret, start, hwSrrTot); + pack_into_bits(ret, start, hwMeanZ); + pack_into_bits(ret, start, hwHoe); return ret; } inline static HadCaloObj unpack(const ap_uint &src) { @@ -59,6 +72,9 @@ namespace l1ct { unpack_from_bits(src, start, ret.hwPhi); unpack_from_bits(src, start, ret.hwEmPt); unpack_from_bits(src, start, ret.hwEmID); + unpack_from_bits(src, start, ret.hwSrrTot); + unpack_from_bits(src, start, ret.hwMeanZ); + unpack_from_bits(src, start, ret.hwHoe); return ret; } }; @@ -70,10 +86,13 @@ namespace l1ct { eta_t hwEta; // relative to the region center, at calo phi_t hwPhi; // relative to the region center, at calo emid_t hwEmID; + srrtot_t hwSrrTot; + meanz_t hwMeanZ; + hoe_t hwHoe; inline bool operator==(const EmCaloObj &other) const { return hwPt == other.hwPt && hwEta == other.hwEta && hwPhi == other.hwPhi && hwPtErr == other.hwPtErr && - hwEmID == other.hwEmID; + hwEmID == other.hwEmID && hwSrrTot == other.hwSrrTot && hwMeanZ == other.hwMeanZ && hwHoe == other.hwHoe; } inline bool operator>(const EmCaloObj &other) const { return hwPt > other.hwPt; } @@ -85,6 +104,9 @@ namespace l1ct { hwEta = 0; hwPhi = 0; hwEmID = 0; + hwSrrTot = 0; + hwMeanZ = 0; + hwHoe = 0; } int intPt() const { return Scales::intPt(hwPt); } @@ -95,8 +117,12 @@ namespace l1ct { float floatPtErr() const { return Scales::floatPt(hwPtErr); } float floatEta() const { return Scales::floatEta(hwEta); } float floatPhi() const { return Scales::floatPhi(hwPhi); } + float floatSrrTot() const { return Scales::floatSrrTot(hwSrrTot); }; + float floatMeanZ() const { return Scales::floatMeanZ(hwMeanZ); }; + float floatHoe() const { return Scales::floatHoe(hwHoe); }; - static const int BITWIDTH = pt_t::width + eta_t::width + phi_t::width + pt_t::width + emid_t::width; + static const int BITWIDTH = pt_t::width + eta_t::width + phi_t::width + pt_t::width + emid_t::width + + srrtot_t::width + meanz_t::width + hoe_t::width; inline ap_uint pack() const { ap_uint ret; unsigned int start = 0; @@ -105,6 +131,9 @@ namespace l1ct { pack_into_bits(ret, start, hwPhi); pack_into_bits(ret, start, hwPtErr); pack_into_bits(ret, start, hwEmID); + pack_into_bits(ret, start, hwSrrTot); + pack_into_bits(ret, start, hwMeanZ); + pack_into_bits(ret, start, hwHoe); return ret; } inline static EmCaloObj unpack(const ap_uint &src) { @@ -115,6 +144,10 @@ namespace l1ct { unpack_from_bits(src, start, ret.hwPhi); unpack_from_bits(src, start, ret.hwPtErr); unpack_from_bits(src, start, ret.hwEmID); + unpack_from_bits(src, start, ret.hwSrrTot); + unpack_from_bits(src, start, ret.hwMeanZ); + unpack_from_bits(src, start, ret.hwHoe); + return ret; } }; @@ -130,6 +163,9 @@ namespace l1ct { z0_t hwZ0; dxy_t hwDxy; tkquality_t hwQuality; + // FIXME: these variables are actually not in the track word or not in this format... + stub_t hwStubs; + chi2_t hwChi2; enum TkQuality { PFLOOSE = 1, PFTIGHT = 2 }; bool isPFLoose() const { return hwQuality[0]; } @@ -140,7 +176,7 @@ namespace l1ct { inline bool operator==(const TkObj &other) const { return hwPt == other.hwPt && hwEta == other.hwEta && hwPhi == other.hwPhi && hwDEta == other.hwDEta && hwDPhi == other.hwDPhi && hwZ0 == other.hwZ0 && hwDxy == other.hwDxy && hwCharge == other.hwCharge && - hwQuality == other.hwQuality; + hwQuality == other.hwQuality && hwChi2 == other.hwChi2 && hwStubs == other.hwStubs; } inline bool operator>(const TkObj &other) const { return hwPt > other.hwPt; } @@ -156,6 +192,8 @@ namespace l1ct { hwDxy = 0; hwCharge = false; hwQuality = 0; + hwStubs = 0; + hwChi2 = 0; } int intPt() const { return Scales::intPt(hwPt); } @@ -173,9 +211,10 @@ namespace l1ct { float floatVtxPhi() const { return Scales::floatPhi(hwVtxPhi()); } float floatZ0() const { return Scales::floatZ0(hwZ0); } float floatDxy() const { return Scales::floatDxy(hwDxy); } + float floatChi2() const { return Scales::floatChi2(hwChi2); } static const int BITWIDTH = pt_t::width + eta_t::width + phi_t::width + tkdeta_t::width + tkdphi_t::width + 1 + - z0_t::width + dxy_t::width + tkquality_t::width; + z0_t::width + dxy_t::width + tkquality_t::width + stub_t::width + chi2_t::width; inline ap_uint pack() const { ap_uint ret; unsigned int start = 0; @@ -188,6 +227,8 @@ namespace l1ct { pack_into_bits(ret, start, hwZ0); pack_into_bits(ret, start, hwDxy); pack_into_bits(ret, start, hwQuality); + pack_into_bits(ret, start, hwStubs); + pack_into_bits(ret, start, hwChi2); return ret; } inline static TkObj unpack(const ap_uint &src) { @@ -202,6 +243,8 @@ namespace l1ct { unpack_from_bits(src, start, ret.hwZ0); unpack_from_bits(src, start, ret.hwDxy); unpack_from_bits(src, start, ret.hwQuality); + unpack_from_bits(src, start, ret.hwStubs); + unpack_from_bits(src, start, ret.hwChi2); return ret; } }; diff --git a/DataFormats/L1TParticleFlow/src/classes_def.xml b/DataFormats/L1TParticleFlow/src/classes_def.xml index ff532e682e559..c0ec668e42d2e 100644 --- a/DataFormats/L1TParticleFlow/src/classes_def.xml +++ b/DataFormats/L1TParticleFlow/src/classes_def.xml @@ -1,6 +1,7 @@ - + + diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/conifer.h b/L1Trigger/Phase2L1ParticleFlow/interface/conifer.h index eb258b6088127..5b53956a4d032 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/conifer.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/conifer.h @@ -1,6 +1,6 @@ #ifndef L1TRIGGER_PHASE2L1PARTICLEFLOW_CONNIFER_H #define L1TRIGGER_PHASE2L1PARTICLEFLOW_CONNIFER_H -#include "FWCore/Utilities/interface/Exception.h" +//#include "FWCore/Utilities/interface/Exception.h" #include "nlohmann/json.hpp" #include @@ -115,16 +115,16 @@ namespace conifer { std::vector decision_function(std::vector x) const { /* Do the prediction */ - if (x.size() != n_features) { - throw cms::Exception("RuntimeError") - << "Conifer : Size of feature vector mismatches expected n_features" << std::endl; - } + //if (x.size() != n_features) { + // throw cms::Exception("RuntimeError") + // << "Conifer : Size of feature vector mismatches expected n_features" << std::endl; + //} std::vector values; std::vector> values_trees; values_trees.resize(n_classes); values.resize(n_classes, U(0)); for (unsigned int i = 0; i < n_classes; i++) { - std::transform(trees.begin(), trees.end(), std::back_inserter(values_trees.at(i)), [&i, &x](auto tree_v) { + std::transform(trees.begin(), trees.end(), std::back_inserter(values_trees.at(i)), [&i, &x](std::vector> tree_v) { return tree_v.at(i).decision_function(x); }); if (useAddTree) { @@ -152,4 +152,4 @@ namespace conifer { } // namespace conifer -#endif \ No newline at end of file +#endif diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h index 93a86c5f4fb8d..863e10de33b64 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h @@ -31,6 +31,8 @@ namespace l1ct { std::vector dEtaValues; std::vector dPhiValues; float trkQualityPtMin; // GeV + bool doCompositeTkEle; + unsigned int nCOMPCAND_PER_CLUSTER; bool writeEgSta; struct IsoParameters { @@ -55,7 +57,7 @@ namespace l1ct { EGIsoEleObjEmu::IsoType hwIsoTypeTkEle; EGIsoObjEmu::IsoType hwIsoTypeTkEm; - bool doCompositeTkEle; + //bool doCompositeTkEle; struct CompIDParameters { CompIDParameters(const edm::ParameterSet &); CompIDParameters(double hoeMin, double hoeMax, double tkptMin, double tkptMax, double srrtotMin, double srrtotMax, double detaMin, double detaMax, double dptMin, double dptMax, double meanzMin, double meanzMax, double dphiMin, double dphiMax, double tkchi2Min, double tkchi2Max, double tkz0Min, double tkz0Max, double tknstubsMin, double tknstubsMax, double BDTcut_wp97p5, double BDTcut_wp95p0) @@ -94,7 +96,7 @@ namespace l1ct { double BDTcut_wp95p0; }; - CompIDParameters myCompIDparams; + CompIDParameters compIDparams; int debug = 0; @@ -115,6 +117,8 @@ namespace l1ct { const std::vector &dEtaValues = {0.015, 0.01}, const std::vector &dPhiValues = {0.07, 0.07}, float trkQualityPtMin = 10., + bool doCompositeTkEle = false, + unsigned int nCompCandPerCluster = 4, bool writeEgSta = false, const IsoParameters &tkIsoParams_tkEle = {2., 0.6, 0.03, 0.2}, const IsoParameters &tkIsoParams_tkEm = {2., 0.6, 0.07, 0.3}, @@ -124,7 +128,7 @@ namespace l1ct { bool doPfIso = false, EGIsoEleObjEmu::IsoType hwIsoTypeTkEle = EGIsoEleObjEmu::IsoType::TkIso, EGIsoObjEmu::IsoType hwIsoTypeTkEm = EGIsoObjEmu::IsoType::TkIsoPV, - bool doCompositeTkEle = false, + // FIXME: maybe we round these? const CompIDParameters &myCompIDparams = {-1.0, 1566.547607421875, 1.9501149654388428, 11102.0048828125, 0.0, 0.01274710614234209, -0.24224889278411865, 0.23079538345336914, 0.010325592942535877, 184.92538452148438, 325.0653991699219, 499.6089782714844, -6.281332015991211, 6.280326843261719, 0.024048099294304848, 1258.37158203125, -14.94140625, 14.94140625, 4.0, 6.0, 0.5406244, 0.9693441}, int debug = 0) @@ -147,6 +151,8 @@ namespace l1ct { dEtaValues(dEtaValues), dPhiValues(dPhiValues), trkQualityPtMin(trkQualityPtMin), + doCompositeTkEle(doCompositeTkEle), + nCOMPCAND_PER_CLUSTER(nCompCandPerCluster), writeEgSta(writeEgSta), tkIsoParams_tkEle(tkIsoParams_tkEle), tkIsoParams_tkEm(tkIsoParams_tkEm), @@ -156,8 +162,7 @@ namespace l1ct { doPfIso(doPfIso), hwIsoTypeTkEle(hwIsoTypeTkEle), hwIsoTypeTkEm(hwIsoTypeTkEm), - doCompositeTkEle(doCompositeTkEle), - myCompIDparams(myCompIDparams), + compIDparams(myCompIDparams), debug(debug) {} }; @@ -383,7 +388,7 @@ namespace l1ct { z0_t z0) const; PFTkEGAlgoEmuConfig cfg; - std::unique_ptr,ap_fixed<22,3,AP_RND_CONV,AP_SAT>,0>> composite_bdt_; + conifer::BDT,ap_fixed<22,3,AP_RND_CONV,AP_SAT>,0> * composite_bdt_; int debug_; }; } // namespace l1ct diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc index 9cade17b5c882..c2083cb1a04e3 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc @@ -38,6 +38,7 @@ #include "DataFormats/L1Trigger/interface/EGamma.h" #include "DataFormats/L1TCorrelator/interface/TkEm.h" #include "DataFormats/L1TCorrelator/interface/TkEmFwd.h" +#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" //-------------------------------------------------------------------------------------------------- class L1TCorrelatorLayer1Producer : public edm::stream::EDProducer<> { @@ -72,7 +73,6 @@ class L1TCorrelatorLayer1Producer : public edm::stream::EDProducer<> { std::unique_ptr l1tkegalgo_; std::unique_ptr l1tkegsorter_; - bool writeEgSta_; // Region dump const std::string regionDumpName_; bool writeRawHgcalCluster_; @@ -176,7 +176,7 @@ L1TCorrelatorLayer1Producer::L1TCorrelatorLayer1Producer(const edm::ParameterSet #if 0 // LATER produces("TKVtx"); #endif -#if 1 // LATER +#if 0 // LATER produces>("DecodedTK"); #endif @@ -369,7 +369,7 @@ void L1TCorrelatorLayer1Producer::produce(edm::Event &iEvent, const edm::EventSe iEvent.put(fetchHadCalo(), "Calo"); iEvent.put(fetchTracks(), "TK"); - #if 1 + #if 0 iEvent.put(fetchDecodedTracks(), "DecodedTK"); #endif @@ -643,7 +643,7 @@ void L1TCorrelatorLayer1Producer::addDecodedTrack(l1ct::DetectorSector 0; } // CMSSW-only extra info - tkAndSel.first.hwChi2 = round(t.chi2() * 10); + tkAndSel.first.hwChi2 = l1ct::Scales::makeChi2(t.chi2()); tkAndSel.first.hwStubs = t.nStubs(); tkAndSel.first.simPt = t.pt(); tkAndSel.first.simCaloEta = t.caloEta(); @@ -653,6 +653,7 @@ void L1TCorrelatorLayer1Producer::addDecodedTrack(l1ct::DetectorSector &sec, const l1t::PFCluster &c) { l1ct::EmCaloObjEmu calo; + // set the endcap-sepcific variables to default value: + calo.clear(); calo.hwPt = l1ct::Scales::makePtFromFloat(c.pt()); calo.hwEta = l1ct::Scales::makeGlbEta(c.eta()) - sec.region.hwEtaCenter; // important to enforce that the region boundary is on a discrete value diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/PFClusterProducerFromHGC3DClusters.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/PFClusterProducerFromHGC3DClusters.cc index 7d393fe347bd9..8594fc19f3c16 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/PFClusterProducerFromHGC3DClusters.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/PFClusterProducerFromHGC3DClusters.cc @@ -156,6 +156,10 @@ void l1tpf::PFClusterProducerFromHGC3DClusters::produce(edm::Event &iEvent, cons corrector_.correctPt(cluster); cluster.setPtError(resol_(cluster.pt(), std::abs(cluster.eta()))); + // We se the cluster shape variables used downstream + cluster.setAbsZBarycenter(fabs(it->zBarycenter())); + cluster.setSigmaRR(it->sigmaRRTot()); + out->push_back(cluster); out->back().addConstituent(edm::Ptr(multiclusters, multiclusters->key(it))); if (hasEmId_) { diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py index 922fc245e81d2..657f51abd5b91 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py @@ -52,6 +52,7 @@ hwIsoTypeTkEle=cms.uint32(0), hwIsoTypeTkEm=cms.uint32(2), doCompositeTkEle=cms.bool(False), + nCOMPCAND_PER_CLUSTER=cms.uint32(3), compositeParametersTkEle=cms.PSet( # Parameters used to normalize input features hoeMin=cms.double(-1.0), hoeMax=cms.double(1566.547607421875), diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py index 686995437e5b3..abcf1200d706f 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py @@ -262,7 +262,7 @@ doEndcapHwQual=True, writeBeforeBremRecovery=False, writeEGSta=True, - doCompositeTkEle=True + doCompositeTkEle=True, trkQualityPtMin=0.), # This should be 10 GeV when doCompositeTkEle=False tkEgSorterParameters=tkEgSorterParameters.clone( nObjToSort = 5 diff --git a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pfeginput_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pfeginput_ref.cpp index 71716e7f542c2..560e0eb96df7f 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pfeginput_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pfeginput_ref.cpp @@ -31,6 +31,10 @@ void EGInputSelectorEmulator::select_eginput(const l1ct::HadCaloObjEmu &in, out.hwPtErr = 0; // shift to get rid of PFEM ID bit (more usable final EG quality) out.hwEmID = (in.hwEmID >> 1); + + out.hwSrrTot = in.hwSrrTot; + out.hwMeanZ = in.hwMeanZ; + out.hwHoe = in.hwHoe; valid_out = (in.hwEmID & cfg.idMask) != 0; } diff --git a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp index ff7f54cfc7135..c0f7f47964b03 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp @@ -7,9 +7,12 @@ #include #include #include +#include + +#include "DataFormats/L1TParticleFlow/interface/PFTrack.h" +#include "DataFormats/L1TParticleFlow/interface/PFCluster.h" + -#include "L1Trigger/Phase2L1ParticleFlow/src/dbgPrintf.h" -#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" using namespace l1ct; #ifdef CMSSW_GIT_HASH @@ -32,6 +35,8 @@ l1ct::PFTkEGAlgoEmuConfig::PFTkEGAlgoEmuConfig(const edm::ParameterSet &pset) dEtaValues(pset.getParameter>("dEtaValues")), dPhiValues(pset.getParameter>("dPhiValues")), trkQualityPtMin(pset.getParameter("trkQualityPtMin")), + doCompositeTkEle(pset.getParameter("doCompositeTkEle")), + nCOMPCAND_PER_CLUSTER(pset.getParameter("nCOMPCAND_PER_CLUSTER")), writeEgSta(pset.getParameter("writeEGSta")), tkIsoParams_tkEle(pset.getParameter("tkIsoParametersTkEle")), tkIsoParams_tkEm(pset.getParameter("tkIsoParametersTkEm")), @@ -41,8 +46,7 @@ l1ct::PFTkEGAlgoEmuConfig::PFTkEGAlgoEmuConfig(const edm::ParameterSet &pset) doPfIso(pset.getParameter("doPfIso")), hwIsoTypeTkEle(static_cast(pset.getParameter("hwIsoTypeTkEle"))), hwIsoTypeTkEm(static_cast(pset.getParameter("hwIsoTypeTkEm"))), - doCompositeTkEle(pset.getParameter("doCompositeTkEle")), - myCompIDparams(pset.getParameter("compositeParametersTkEle")), + compIDparams(pset.getParameter("compositeParametersTkEle")), debug(pset.getUntrackedParameter("debug", 0)) {} l1ct::PFTkEGAlgoEmuConfig::IsoParameters::IsoParameters(const edm::ParameterSet &pset) @@ -82,8 +86,14 @@ composite_bdt_(nullptr), debug_(cfg.debug) { if(cfg.doCompositeTkEle) { //FIXME: make the name of the file configurable - auto resolvedFileName = edm::FileInPath("L1Trigger/Phase2L1ParticleFlow/data/compositeID.json").fullPath(); - composite_bdt_ = std::make_unique,ap_fixed<22,3,AP_RND_CONV,AP_SAT>,0>> (resolvedFileName); +#ifdef CMSSW_GIT_HASH + auto resolvedFileName = edm::FileInPath("L1Trigger/Phase2L1ParticleFlow/data/compositeID.json").fullPath(); +#else + auto resolvedFileName = "compositeID.json"; +#endif + std::cout<,ap_fixed<22,3,AP_RND_CONV,AP_SAT>,0> (resolvedFileName); + std::cout<<"declared bdt"< &track, std::vector &emCalo2tk, std::vector &emCaloTkBdtScore) const { - //FIXME: should be configurable - const int nCAND_PER_CLUSTER = 4; unsigned int nTrackMax = std::min(track.size(), cfg.nTRACK_EGIN); + std::cout<<"doing loose dR matching"< candidates; for (unsigned int itk = 0; itk < nTrackMax; ++itk) { + std::cout<<"track "< bool { return a.dpt < b.dpt; }); - unsigned int nCandPerCluster = std::min(candidates.size(), nCAND_PER_CLUSTER); - std::cout << "# composit candidates: " << nCandPerCluster << std::endl; + unsigned int nCandPerCluster = std::min(candidates.size(), cfg.nCOMPCAND_PER_CLUSTER); + std::cout << "# composite candidates: " << nCandPerCluster << std::endl; if(nCandPerCluster == 0) continue; - float bdtWP_MVA = cfg.myCompIDparams.BDTcut_wp97p5; + float bdtWP_MVA = cfg.compIDparams.BDTcut_wp97p5; float bdtWP_XGB = 1. / (1. + std::sqrt((1. - bdtWP_MVA) / (1. + bdtWP_MVA))); // Convert WP value from ROOT.TMVA to XGboost float maxScore = -999; int ibest = -1; for(unsigned int icand = 0; icand < nCandPerCluster; icand++) { auto &cand = candidates[icand]; std::vector emcalo_sel = emcalo; - float score = compute_composite_score(cand, emcalo_sel, track, cfg.myCompIDparams); + float score = compute_composite_score(cand, emcalo_sel, track, cfg.compIDparams); if(score > maxScore) { // if((score > bdtWP_XGB) && (score > maxScore)) { maxScore = score; @@ -255,37 +267,23 @@ float PFTkEGAlgoEmulator::compute_composite_score(CompositeCandidate &cand, const auto &calo = emcalo[cand.cluster_idx]; const auto &tk = track[cand.track_idx]; - // FIXME: using these two floats crashes code... - float srrtot_f = dynamic_cast(calo.src->constituentsAndFractions()[0].first.get())->sigmaRRTot(); - float meanz_f = abs(dynamic_cast(calo.src->constituentsAndFractions()[0].first.get())->zBarycenter()); - // Call and normalize input feature values, then cast to ap_fixed. // Note that for some features (e.g. track pT) we call the floating point representation, but that's already quantized! // Several other features, such as chi2 or most cluster features, are not quantized before casting them to ap_fixed. - ap_fixed<22,3,AP_RND_CONV,AP_SAT> hoe = (calo.src->hOverE()-params.hoeMin)/(params.hoeMax-params.hoeMin); + ap_fixed<22,3,AP_RND_CONV,AP_SAT> hoe = (calo.floatHoe()-params.hoeMin)/(params.hoeMax-params.hoeMin); ap_fixed<22,3,AP_RND_CONV,AP_SAT> tkpt = (tk.floatPt()-params.tkptMin)/(params.tkptMax-params.tkptMin); - ap_fixed<22,3,AP_RND_CONV,AP_SAT> srrtot = (srrtot_f-params.srrtotMin)/(params.srrtotMax-params.srrtotMin); + ap_fixed<22,3,AP_RND_CONV,AP_SAT> srrtot = (calo.floatSrrTot()-params.srrtotMin)/(params.srrtotMax-params.srrtotMin); ap_fixed<22,3,AP_RND_CONV,AP_SAT> deta = (tk.floatEta() - calo.floatEta()-params.detaMin)/(params.detaMax-params.detaMin); // FIXME: do we really need dpt to be a ratio? ap_fixed<22,3,AP_RND_CONV,AP_SAT> dpt = ((tk.floatPt()/calo.floatPt())-params.dptMin)/(params.dptMax-params.dptMin); - ap_fixed<22,3,AP_RND_CONV,AP_SAT> meanz = (meanz_f-params.meanzMin)/(params.meanzMax-params.meanzMin); + ap_fixed<22,3,AP_RND_CONV,AP_SAT> meanz = (calo.floatMeanZ()-params.meanzMin)/(params.meanzMax-params.meanzMin); ap_fixed<22,3,AP_RND_CONV,AP_SAT> dphi = (deltaPhi(tk.floatPhi(), calo.floatPhi()) -params.dphiMin)/(params.dphiMax-params.dphiMin); - ap_fixed<22,3,AP_RND_CONV,AP_SAT> chi2 = (tk.src->chi2()-params.tkchi2Min)/(params.tkchi2Max-params.tkchi2Min); + ap_fixed<22,3,AP_RND_CONV,AP_SAT> chi2 = (tk.floatChi2()-params.tkchi2Min)/(params.tkchi2Max-params.tkchi2Min); ap_fixed<22,3,AP_RND_CONV,AP_SAT> tkz0 = (tk.floatZ0()-params.tkz0Min)/(params.tkz0Max-params.tkz0Min); - ap_fixed<22,3,AP_RND_CONV,AP_SAT> nstubs = (tk.src->nStubs()-params.tknstubsMin)/(params.tknstubsMax-params.tknstubsMin); - // std::cout<<"hoe\t"<hOverE()<<"\t"<<(calo.src->hOverE()-params.hoeMin)/(params.hoeMax-params.hoeMin)<chi2()<<"\t"<<(tk.src->chi2()-params.tkchi2Min)/(params.tkchi2Max-params.tkchi2Min)<nStubs()<<"\t"<<(tk.src->nStubs()-params.tknstubsMin)/(params.tknstubsMax-params.tknstubsMin)< nstubs = (tk.hwStubs-params.tknstubsMin)/(params.tknstubsMax-params.tknstubsMin); + // Run BDT inference - vector> inputs = { hoe, tkpt, srrtot, deta, dpt, meanz, dphi, chi2, tkz0, nstubs } ; + std::vector> inputs = { hoe, tkpt, srrtot, deta, dpt, meanz, dphi, chi2, tkz0, nstubs } ; auto bdt_score = composite_bdt_->decision_function(inputs); float bdt_score_CON = bdt_score[0]; @@ -321,7 +319,7 @@ void PFTkEGAlgoEmulator::run(const PFInputRegion &in, OutputRegion &out) const { << std::endl; } } - + std::cout<<"running"< emCalo2tk(emcalo_sel.size(), -1); std::vector emCaloTkBdtScore(emcalo_sel.size(), -999); + std::cout<<"about to start matching"< out.hwEmPt = w_empt * l1ct::pt_t(l1ct::Scales::INTPT_LSB); out.hwEmID = w_qual; + // FIXME: cluster-shape variables use by composite-ID need to be added here + return out; } diff --git a/L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/compositeID.json b/L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/compositeID.json deleted file mode 100644 index c7cfb823abdd1..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/compositeID.json +++ /dev/null @@ -1 +0,0 @@ -{"max_depth": 3, "n_trees": 10, "n_classes": 2, "n_features": 10, "trees": [[{"feature": [0, 2, 0, 1, 1, 1, 0, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000703328056, 0.332558006, 0.000770003942, 0.000253851322, 0.000915614248, 0.000717727817, 0.000839613494, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.0352521278, 0.590574145, -0.54147464, 0.185335979, -0.503847778, 0.257735074, -0.503176689, -0.595144987], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [0, 2, 0, 1, 1, 1, 5, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000697391108, 0.340182066, 0.000747981714, 0.000371906266, 0.00168169127, 0.000789736281, 0.0738627762, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.0733107328, 0.458879024, -0.396272093, 0.309196174, -0.360979736, 0.275205106, -0.349665552, -0.460340261], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [0, 1, 1, 4, 2, 5, 5, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000694462447, 0.000501940958, 0.00114817475, 0.000432471978, 0.370490968, 0.0738627762, 0.0940087885, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.380878597, -0.283362597, 0.397695929, -0.294095337, -0.280492604, -0.399287552, 0.310794801, -0.362855673], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [0, 2, 0, 3, 1, 2, 0, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000709932996, 0.321239412, 0.000803797971, 0.602574646, 0.00082699745, 0.336098373, 0.000910258619, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.353067011, -0.383320004, -0.345221847, 0.114074722, 0.156569079, -0.339963675, -0.309075147, -0.367914289], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [0, 1, 1, 4, 2, 7, 2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000687893422, 0.000643459614, 0.000599945546, 0.000337057747, 0.358681321, 0.0250910092, 0.343107432, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.349036545, -0.16346851, 0.341645479, -0.123347722, -0.342872322, -0.239095137, 0.228956148, -0.307008922], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [0, 1, 0, 4, 2, 1, 0, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000722329598, 0.000350446324, 0.000839613494, 0.000351717492, 0.365002632, 0.00043548242, 0.000917885685, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.253911734, -0.244474351, 0.316286922, -0.238444999, -0.279888034, 0.0941597223, -0.259838462, -0.331297874], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [0, 3, 5, 8, 3, 1, 1, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000677733915, 0.438798308, 0.0826261193, 0.850980401, 0.592054725, 0.00032415273, 0.00283236895, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, -0.413880289, 0.257492244, 0.303176314, -0.260525465, -0.260994583, 0.136501729, -0.314218611, 0.153914765], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [1, 0, 2, 2, 0, 5, 2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.00122120825, 0.000803797971, 0.33025986, 0.301885217, 0.000917885685, 0.133645415, 0.387662947, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.145611659, -0.207396939, -0.205260798, -0.312455952, 0.311472386, -0.278300524, -0.0131324222, -0.31602335], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [1, 5, 2, 2, 5, 0, 4, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.00148988748, 0.0920395404, 0.357989103, 0.310617983, 0.116608188, 0.00109421043, 0.00491232192, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.13668026, -0.192945674, -0.225791544, -0.311047494, 0.302013129, -0.30914408, -0.0594445392, -0.352083832], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}], [{"feature": [1, 7, 2, 4, 0, 5, 2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.000719727308, 0.0146044344, 0.315754294, 0.000207430043, 0.000914971635, 0.114892721, 0.404698849, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.209280849, -0.279667407, 0.159891739, -0.288006365, 0.288983703, -0.244072288, 0.0150486836, -0.272770345], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3], "iLeaf": [7, 8, 9, 10, 11, 12, 13, 14]}]], "init_predict": [0, 0], "norm": 1} \ No newline at end of file diff --git a/L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/conifer.h b/L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/conifer.h deleted file mode 100644 index d9529aa9dcb90..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/conifer.h +++ /dev/null @@ -1,158 +0,0 @@ -#ifndef CONIFER_CPP_H__ -#define CONIFER_CPP_H__ -#include "nlohmann/json.hpp" -#include -#include - -namespace conifer{ - -/* --- -* Balanced tree reduce implementation. -* Reduces an array of inputs to a single value using the template binary operator 'Op', -* for example summing all elements with Op_add, or finding the maximum with Op_max -* Use only when the input array is fully unrolled. Or, slice out a fully unrolled section -* before applying and accumulate the result over the rolled dimension. -* Required for emulation to guarantee equality of ordering. -* --- */ -constexpr int floorlog2(int x) { return (x < 2) ? 0 : 1 + floorlog2(x / 2); } - -template -constexpr int pow(int x) { - return x == 0 ? 1 : B * pow(x - 1); -} - -constexpr int pow2(int x) { return pow<2>(x); } - -template -T reduce(std::vector x, Op op) { - int N = x.size(); - int leftN = pow2(floorlog2(N - 1)) > 0 ? pow2(floorlog2(N - 1)) : 0; - //static constexpr int rightN = N - leftN > 0 ? N - leftN : 0; - if (N == 1) { - return x.at(0); - } else if (N == 2) { - return op(x.at(0), x.at(1)); - } else { - std::vector left(x.begin(), x.begin() + leftN); - std::vector right(x.begin() + leftN, x.end()); - return op(reduce(left, op), reduce(right, op)); - } -} - -template -class OpAdd { -public: - T operator()(T a, T b) { return a + b; } -}; - -template -class DecisionTree{ - -private: - std::vector feature; - std::vector children_left; - std::vector children_right; - std::vector threshold_; - std::vector value_; - std::vector threshold; - std::vector value; - -public: - - U decision_function(std::vector x) const{ - /* Do the prediction */ - int i = 0; - while(feature[i] != -2){ // continue until reaching leaf - bool comparison = x[feature[i]] <= threshold_[i]; - i = comparison ? children_left[i] : children_right[i]; - } - return value_[i]; - } - - void init_(){ - /* Since T, U types may not be readable from the JSON, read them to double and the cast them here */ - std::transform(threshold.begin(), threshold.end(), std::back_inserter(threshold_), - [](double t) -> T { return (T) t; }); - std::transform(value.begin(), value.end(), std::back_inserter(value_), - [](double v) -> U { return (U) v; }); - } - - // Define how to read this class to/from JSON - NLOHMANN_DEFINE_TYPE_INTRUSIVE(DecisionTree, feature, children_left, children_right, threshold, value); - -}; // class DecisionTree - -template -class BDT{ - -private: - - int n_classes; - int n_trees; - unsigned int n_features; - std::vector init_predict; - std::vector init_predict_; - // vector of decision trees: outer dimension tree, inner dimension class - std::vector>> trees; - OpAdd add; - -public: - - // Define how to read this class to/from JSON - NLOHMANN_DEFINE_TYPE_INTRUSIVE(BDT, n_classes, n_trees, n_features, init_predict, trees); - - BDT(std::string filename){ - /* Construct the BDT from conifer cpp backend JSON file */ - std::ifstream ifs(filename); - nlohmann::json j = nlohmann::json::parse(ifs); - from_json(j, *this); - /* Do some transformation to initialise things into the proper emulation T, U types */ - if(n_classes == 2) n_classes = 1; - std::transform(init_predict.begin(), init_predict.end(), std::back_inserter(init_predict_), - [](double ip) -> U { return (U) ip; }); - for(int i = 0; i < n_trees; i++){ - for(int j = 0; j < n_classes; j++){ - trees.at(i).at(j).init_(); - } - } - } - - std::vector decision_function(std::vector x) const{ - /* Do the prediction */ - // std::cout< values; - std::vector> values_trees; - values_trees.resize(n_classes); - values.resize(n_classes, U(0)); - for(int i = 0; i < n_classes; i++){ - std::transform(trees.begin(), trees.end(), std::back_inserter(values_trees.at(i)), - [&i, &x](auto tree_v){ return tree_v.at(i).decision_function(x); }); - if(useAddTree){ - values.at(i) = init_predict_.at(i); - values.at(i) += reduce>(values_trees.at(i), add); - }else{ - values.at(i) = std::accumulate(values_trees.at(i).begin(), values_trees.at(i).end(), U(init_predict_.at(i))); - } - } - - return values; - } - - std::vector _decision_function_double(std::vector x) const{ - /* Do the prediction with data in/out as double, cast to T, U before prediction */ - std::vector xt; - std::transform(x.begin(), x.end(), std::back_inserter(xt), - [](double xi) -> T { return (T) xi; }); - std::vector y = decision_function(xt); - std::vector yd; - std::transform(y.begin(), y.end(), std::back_inserter(yd), - [](U yi) -> double { return (double) yi; }); - return yd; - } - -}; // class BDT - -} // namespace conifer - -#endif \ No newline at end of file From 03daa243ba152da1118c5931937c0a5ff252c34a Mon Sep 17 00:00:00 2001 From: Peter Eduard Meiring Date: Thu, 10 Nov 2022 19:08:09 +0100 Subject: [PATCH 04/13] change to new model with max tree depth of 4 --- .../Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h | 2 +- .../Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py | 4 ++-- L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h index 863e10de33b64..16b4aad96bede 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h @@ -129,7 +129,7 @@ namespace l1ct { EGIsoEleObjEmu::IsoType hwIsoTypeTkEle = EGIsoEleObjEmu::IsoType::TkIso, EGIsoObjEmu::IsoType hwIsoTypeTkEm = EGIsoObjEmu::IsoType::TkIsoPV, // FIXME: maybe we round these? - const CompIDParameters &myCompIDparams = {-1.0, 1566.547607421875, 1.9501149654388428, 11102.0048828125, 0.0, 0.01274710614234209, -0.24224889278411865, 0.23079538345336914, 0.010325592942535877, 184.92538452148438, 325.0653991699219, 499.6089782714844, -6.281332015991211, 6.280326843261719, 0.024048099294304848, 1258.37158203125, -14.94140625, 14.94140625, 4.0, 6.0, 0.5406244, 0.9693441}, + const CompIDParameters &myCompIDparams = {-1.0, 1566.547607421875, 1.9501149654388428, 11102.0048828125, 0.0, 0.01274710614234209, -0.24224889278411865, 0.23079538345336914, 0.010325592942535877, 184.92538452148438, 325.0653991699219, 499.6089782714844, -6.281332015991211, 6.280326843261719, 0.024048099294304848, 1258.37158203125, -14.94140625, 14.94140625, 4.0, 6.0, 0.7927004, 0.9826955}, int debug = 0) : nTRACK(nTrack), diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py index 657f51abd5b91..e270bb99bca73 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py @@ -74,8 +74,8 @@ tkz0Max=cms.double(14.94140625), tknstubsMin=cms.double(4.0), tknstubsMax=cms.double(6.0), - BDTcut_wp97p5=cms.double(0.5406244), - BDTcut_wp95p0=cms.double(0.9693441), + BDTcut_wp97p5=cms.double(0.7927004), + BDTcut_wp95p0=cms.double(0.9826955), ), ) diff --git a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp index c0f7f47964b03..d8966b2e9db0c 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp @@ -89,7 +89,7 @@ debug_(cfg.debug) { #ifdef CMSSW_GIT_HASH auto resolvedFileName = edm::FileInPath("L1Trigger/Phase2L1ParticleFlow/data/compositeID.json").fullPath(); #else - auto resolvedFileName = "compositeID.json"; + auto resolvedFileName = "compositeID.json"; #endif std::cout<,ap_fixed<22,3,AP_RND_CONV,AP_SAT>,0> (resolvedFileName); From a08f0ea7cd7afabe6e9cdfd3e056c22da087c1b6 Mon Sep 17 00:00:00 2001 From: CI bot Date: Wed, 23 Nov 2022 17:01:14 +0100 Subject: [PATCH 05/13] add variales for composite ID to HGCAL raw dataformat --- DataFormats/L1TParticleFlow/interface/datatypes.h | 3 ++- .../plugins/L1TCorrelatorLayer1Producer.cc | 14 +++++++++++++- .../src/l1-converters/hgcalinputt_ref.cpp | 10 +++++++++- 3 files changed, 24 insertions(+), 3 deletions(-) diff --git a/DataFormats/L1TParticleFlow/interface/datatypes.h b/DataFormats/L1TParticleFlow/interface/datatypes.h index 719ab378d3205..e579cd8493163 100644 --- a/DataFormats/L1TParticleFlow/interface/datatypes.h +++ b/DataFormats/L1TParticleFlow/interface/datatypes.h @@ -40,7 +40,6 @@ namespace l1ct { typedef ap_uint<13> tk2calo_dq_t; typedef ap_uint<4> egquality_t; typedef ap_uint<3> stub_t; - // FIXME: random choice for the various bitwidth below!!!! typedef ap_ufixed<12, 10, AP_TRN, AP_SAT> chi2_t; typedef ap_ufixed<10, 1, AP_TRN, AP_SAT> srrtot_t; typedef ap_uint<8> meanz_t; // mean - SCALE_MEANZ = 320 @@ -157,6 +156,8 @@ namespace l1ct { constexpr float DXY_LSB = 0.05; constexpr float PUPPIW_LSB = 1.0 / 256; constexpr float MEANZ_SCALE = 320.; + constexpr float SRRTOT_LSB = pow(2,-9); + constexpr float HOE_LSB = pow(2,-5); inline float floatPt(pt_t pt) { return pt.to_float(); } inline float floatPt(dpt_t pt) { return pt.to_float(); } diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc index c2083cb1a04e3..3031ffd06460f 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc @@ -705,13 +705,25 @@ void L1TCorrelatorLayer1Producer::addRawHgcalCluster(l1ct::DetectorSector w_phi = round(sec.region.localPhi(c.phi()) / ETAPHI_LSB); ap_uint<10> w_qual = c.hwQual(); + ap_uint<13> w_srrtot = round(c.sigmaRR()/l1ct::Scales::SRRTOT_LSB); + ap_uint<12> w_meanz = round(c.absZBarycenter()); + ap_uint<12> w_hoe = round(c.hOverE()/l1ct::Scales::HOE_LSB); + cwrd(13, 0) = w_pt; cwrd(27, 14) = w_empt; cwrd(72, 64) = w_eta; cwrd(81, 73) = w_phi; cwrd(115, 106) = w_qual; - // FIXME: cluster-shape variables use by composite-ID need to be added here + // FIXME: we add the variables use by composite-ID. The definitin will have to be reviewd once the + // hgc format is better defined. For now we use + // hwMeanZ = word 1 bits 30-19 + // hwSrrTot = word 3 bits 21 - 9 + // hoe = word 1 bits 63-52 (currently spare) + cwrd(213, 201) = w_srrtot; + cwrd(94, 83) = w_meanz; + // FIXME: we use a spare space in the word for hoe which is not in the current interface + cwrd(127, 116) = w_hoe; sec.obj.push_back(cwrd); } diff --git a/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp index 3a079ed5bd6c4..48d8fc85fb2cc 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp @@ -8,7 +8,12 @@ l1ct::HadCaloObjEmu l1ct::HgcalClusterDecoderEmulator::decode(const ap_uint<256> ap_int<9> w_eta = in(72, 64); ap_int<9> w_phi = in(81, 73); ap_uint<10> w_qual = in(115, 106); + ap_uint<13> w_srrtot = in(213, 201); + ap_uint<12> w_meanz = in(94, 83); + // FIXME: we use a spare space in the word for hoe which is not in the current interface + ap_uint<12> w_hoe = in(127, 116); + l1ct::HadCaloObjEmu out; out.clear(); out.hwPt = w_pt * l1ct::pt_t(l1ct::Scales::INTPT_LSB); @@ -17,7 +22,10 @@ l1ct::HadCaloObjEmu l1ct::HgcalClusterDecoderEmulator::decode(const ap_uint<256> out.hwEmPt = w_empt * l1ct::pt_t(l1ct::Scales::INTPT_LSB); out.hwEmID = w_qual; - // FIXME: cluster-shape variables use by composite-ID need to be added here + // FIXME: variables use by composite-ID need to be added here + out.hwSrrTot = w_srrtot * l1ct::srrtot_t(l1ct::Scales::SRRTOT_LSB); + out.hwMeanZ = l1ct::meanz_t(w_meanz - l1ct::Scales::MEANZ_SCALE); + calo.hwHoe = w_hoe * l1ct::hoe_t(l1ct::Scales::HOE_LSB); return out; } From 4d63273217f4a2db15617810281a6c39224d7195 Mon Sep 17 00:00:00 2001 From: CI bot Date: Mon, 21 Nov 2022 16:38:32 +0100 Subject: [PATCH 06/13] Replace egamma BDT with model using no input scaling --- .../interface/egamma/pftkegalgo_ref.h | 2 +- .../src/egamma/pftkegalgo_ref.cpp | 26 +++++++++---------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h index 16b4aad96bede..65dbf5f0fd864 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h @@ -388,7 +388,7 @@ namespace l1ct { z0_t z0) const; PFTkEGAlgoEmuConfig cfg; - conifer::BDT,ap_fixed<22,3,AP_RND_CONV,AP_SAT>,0> * composite_bdt_; + conifer::BDT,ap_fixed<12,3,AP_RND_CONV,AP_SAT>,0> * composite_bdt_; int debug_; }; } // namespace l1ct diff --git a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp index d8966b2e9db0c..d89f3f989c3e8 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp @@ -92,7 +92,7 @@ debug_(cfg.debug) { auto resolvedFileName = "compositeID.json"; #endif std::cout<,ap_fixed<22,3,AP_RND_CONV,AP_SAT>,0> (resolvedFileName); + composite_bdt_ = new conifer::BDT,ap_fixed<12,3,AP_RND_CONV,AP_SAT>,0> (resolvedFileName); std::cout<<"declared bdt"< hoe = (calo.floatHoe()-params.hoeMin)/(params.hoeMax-params.hoeMin); - ap_fixed<22,3,AP_RND_CONV,AP_SAT> tkpt = (tk.floatPt()-params.tkptMin)/(params.tkptMax-params.tkptMin); - ap_fixed<22,3,AP_RND_CONV,AP_SAT> srrtot = (calo.floatSrrTot()-params.srrtotMin)/(params.srrtotMax-params.srrtotMin); - ap_fixed<22,3,AP_RND_CONV,AP_SAT> deta = (tk.floatEta() - calo.floatEta()-params.detaMin)/(params.detaMax-params.detaMin); + ap_fixed<21,12,AP_RND_CONV,AP_SAT> hoe = calo.floatHoe(); + ap_fixed<21,12,AP_RND_CONV,AP_SAT> tkpt = tk.floatPt(); + ap_fixed<21,12,AP_RND_CONV,AP_SAT> srrtot = calo.floatSrrTot(); + ap_fixed<21,12,AP_RND_CONV,AP_SAT> deta = tk.floatEta() - calo.floatEta(); // FIXME: do we really need dpt to be a ratio? - ap_fixed<22,3,AP_RND_CONV,AP_SAT> dpt = ((tk.floatPt()/calo.floatPt())-params.dptMin)/(params.dptMax-params.dptMin); - ap_fixed<22,3,AP_RND_CONV,AP_SAT> meanz = (calo.floatMeanZ()-params.meanzMin)/(params.meanzMax-params.meanzMin); - ap_fixed<22,3,AP_RND_CONV,AP_SAT> dphi = (deltaPhi(tk.floatPhi(), calo.floatPhi()) -params.dphiMin)/(params.dphiMax-params.dphiMin); - ap_fixed<22,3,AP_RND_CONV,AP_SAT> chi2 = (tk.floatChi2()-params.tkchi2Min)/(params.tkchi2Max-params.tkchi2Min); - ap_fixed<22,3,AP_RND_CONV,AP_SAT> tkz0 = (tk.floatZ0()-params.tkz0Min)/(params.tkz0Max-params.tkz0Min); - ap_fixed<22,3,AP_RND_CONV,AP_SAT> nstubs = (tk.hwStubs-params.tknstubsMin)/(params.tknstubsMax-params.tknstubsMin); + ap_fixed<21,12,AP_RND_CONV,AP_SAT> dpt = (tk.floatPt()/calo.floatPt()); + ap_fixed<21,12,AP_RND_CONV,AP_SAT> meanz = calo.floatMeanZ(); + ap_fixed<21,12,AP_RND_CONV,AP_SAT> dphi = deltaPhi(tk.floatPhi(), calo.floatPhi()); + ap_fixed<21,12,AP_RND_CONV,AP_SAT> chi2 = tk.floatChi2(); + ap_fixed<21,12,AP_RND_CONV,AP_SAT> tkz0 = tk.floatZ0(); + ap_fixed<21,12,AP_RND_CONV,AP_SAT> nstubs = tk.hwStubs; // Run BDT inference - std::vector> inputs = { hoe, tkpt, srrtot, deta, dpt, meanz, dphi, chi2, tkz0, nstubs } ; - auto bdt_score = composite_bdt_->decision_function(inputs); + std::vector> inputs = { hoe, tkpt, srrtot, deta, dpt, meanz, dphi, chi2, tkz0, nstubs } ; + std::vector> bdt_score = composite_bdt_->decision_function(inputs); float bdt_score_CON = bdt_score[0]; float bdt_score_XGB = 1/(1+exp(-bdt_score_CON)); // Map Conifer score to XGboost score. (same as scipy.expit) From 85836bfcd030f1b77beae688c3216264007bbc8a Mon Sep 17 00:00:00 2001 From: CI bot Date: Tue, 22 Nov 2022 14:39:02 +0100 Subject: [PATCH 07/13] Update egamma BDT to model trained on hw values. Use inversion fn from jets --- .../interface/common/inversion.h | 61 +++++++++++++++++++ .../interface/egamma/pftkegalgo_ref.h | 1 + .../src/egamma/pftkegalgo_ref.cpp | 20 +++--- 3 files changed, 72 insertions(+), 10 deletions(-) create mode 100644 L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h b/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h new file mode 100644 index 0000000000000..58a8492bd4def --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h @@ -0,0 +1,61 @@ +#ifndef CC_INVERSION_H__ +#define CC_INVERSION_H__ + +constexpr int ceillog2(int x){ + return (x <= 2) ? 1 : 1 + ceillog2((x+1) / 2); +} + +template +inline float real_val_from_idx(unsigned i){ + // Treat the index as the top N bits + static constexpr int NB = ceillog2(N); // number of address bits for table + data_T x(0); + // The MSB of 1 is implicit in the table + x[x.width-1] = 1; + // So we can use the next NB bits for real data + x(x.width-2, x.width-NB-1) = i; + return (float) x; +} + +template +inline unsigned idx_from_real_val(data_T x){ + // Slice the top N bits to get an index into the table + static constexpr int NB = ceillog2(N); // number of address bits for table + // Slice the top-1 NB bits of the value + // the MSB of '1' is implicit, so only slice below that + ap_uint y = x(x.width-2, x.width-NB-1); + return (unsigned) y(NB-1, 0); +} + +template +void init_invert_table(table_T table_out[N]){ + // The template data_T is the data type used to address the table + for(unsigned i = 0; i < N; i++){ + float x = real_val_from_idx(i); + table_T inv_x = 1 / x; + table_out[i] = inv_x; + } +} + +template +table_t invert_with_shift(in_t in){ + table_t inv_table[N]; + init_invert_table(inv_table); + + // find the first '1' in the denominator + int msb = 0; + for(int b = 0; b < in.width; b++){ + #pragma HLS unroll + if(in[b]) msb = b; + } + // shift up the denominator such that the left-most bit (msb) is '1' + in_t in_shifted = in << (in.width-msb-1); + // lookup the inverse of the shifted input + int idx = idx_from_real_val(in_shifted); + table_t inv_in = inv_table[idx]; + // shift the output back + table_t out = inv_in << (in.width-msb-1); + return out; +} + +#endif \ No newline at end of file diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h index 65dbf5f0fd864..7eb6cf5e449f8 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h @@ -5,6 +5,7 @@ #include "DataFormats/L1TParticleFlow/interface/egamma.h" #include "DataFormats/L1TParticleFlow/interface/pf.h" #include "L1Trigger/Phase2L1ParticleFlow/interface/conifer.h" +#include "L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h" namespace edm { class ParameterSet; diff --git a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp index d89f3f989c3e8..5786ac32bebe0 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp @@ -270,16 +270,16 @@ float PFTkEGAlgoEmulator::compute_composite_score(CompositeCandidate &cand, // Call and normalize input feature values, then cast to ap_fixed. // Note that for some features (e.g. track pT) we call the floating point representation, but that's already quantized! // Several other features, such as chi2 or most cluster features, are not quantized before casting them to ap_fixed. - ap_fixed<21,12,AP_RND_CONV,AP_SAT> hoe = calo.floatHoe(); - ap_fixed<21,12,AP_RND_CONV,AP_SAT> tkpt = tk.floatPt(); - ap_fixed<21,12,AP_RND_CONV,AP_SAT> srrtot = calo.floatSrrTot(); - ap_fixed<21,12,AP_RND_CONV,AP_SAT> deta = tk.floatEta() - calo.floatEta(); - // FIXME: do we really need dpt to be a ratio? - ap_fixed<21,12,AP_RND_CONV,AP_SAT> dpt = (tk.floatPt()/calo.floatPt()); - ap_fixed<21,12,AP_RND_CONV,AP_SAT> meanz = calo.floatMeanZ(); - ap_fixed<21,12,AP_RND_CONV,AP_SAT> dphi = deltaPhi(tk.floatPhi(), calo.floatPhi()); - ap_fixed<21,12,AP_RND_CONV,AP_SAT> chi2 = tk.floatChi2(); - ap_fixed<21,12,AP_RND_CONV,AP_SAT> tkz0 = tk.floatZ0(); + ap_fixed<21,12,AP_RND_CONV,AP_SAT> hoe = calo.hwHoe; + ap_fixed<21,12,AP_RND_CONV,AP_SAT> tkpt = tk.hwPt; + ap_fixed<21,12,AP_RND_CONV,AP_SAT> srrtot = calo.hwSrrTot; + ap_fixed<21,12,AP_RND_CONV,AP_SAT> deta = tk.hwEta - calo.hwEta; + ap_fixed<18,9> calo_invPt = invert_with_shift, 1024>(calo.hwPt); // TODO: this is a guess + ap_fixed<21,12,AP_RND_CONV,AP_SAT> dpt = tk.hwPt * calo_invPt; + ap_fixed<21,12,AP_RND_CONV,AP_SAT> meanz = calo.hwMeanZ; + ap_fixed<21,12,AP_RND_CONV,AP_SAT> dphi = tk.hwPhi - calo.hwPhi; + ap_fixed<21,12,AP_RND_CONV,AP_SAT> chi2 = tk.hwChi2; + ap_fixed<21,12,AP_RND_CONV,AP_SAT> tkz0 = tk.hwZ0; ap_fixed<21,12,AP_RND_CONV,AP_SAT> nstubs = tk.hwStubs; // Run BDT inference From 40a6c93dec977979bb53e8b62fe39ffec56df130 Mon Sep 17 00:00:00 2001 From: Gianluca Date: Wed, 23 Nov 2022 17:15:11 +0100 Subject: [PATCH 08/13] Introduce Slim version of packed objects and handle BDT WP for setting hwQual. --- .../L1TParticleFlow/interface/PFCandidate.h | 1 - .../L1TParticleFlow/interface/bit_encoding.h | 30 ++++ .../L1TParticleFlow/interface/datatypes.h | 20 +-- .../interface/layer1_emulator.h | 5 +- .../L1TParticleFlow/interface/layer1_objs.h | 48 ++++-- .../L1TParticleFlow/src/layer1_emulator.cpp | 12 +- .../HGCalTriggerNtupleHGCMulticlusters.cc | 12 +- .../interface/common/inversion.h | 93 ++++++----- .../Phase2L1ParticleFlow/interface/conifer.h | 27 ++- .../interface/egamma/pftkegalgo_ref.h | 72 +++----- .../interface/l1-converters/tkinput_ref.h | 28 +++- .../plugins/L1TCorrelatorLayer1Producer.cc | 47 +++--- .../PFClusterProducerFromHGC3DClusters.cc | 3 +- .../python/l1TkEgAlgoEmulator_cfi.py | 32 +--- .../python/l1TkEgAlgo_cfi.py | 47 ------ .../python/l1ctLayer1_cff.py | 2 + .../src/egamma/pftkegalgo_ref.cpp | 154 +++++++----------- .../src/l1-converters/hgcalinputt_ref.cpp | 7 +- .../src/l1-converters/tkinput_ref.cpp | 18 +- .../test/make_l1ctLayer1_dumpFiles_cfg.py | 2 +- 20 files changed, 309 insertions(+), 351 deletions(-) delete mode 100644 L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgo_cfi.py diff --git a/DataFormats/L1TParticleFlow/interface/PFCandidate.h b/DataFormats/L1TParticleFlow/interface/PFCandidate.h index 50cf529178dc9..567eb6f80b7fb 100644 --- a/DataFormats/L1TParticleFlow/interface/PFCandidate.h +++ b/DataFormats/L1TParticleFlow/interface/PFCandidate.h @@ -51,7 +51,6 @@ namespace l1t { void setCaloEta(float caloeta) { caloEta_ = caloeta; } void setCaloPhi(float calophi) { caloPhi_ = calophi; } - float z0() const { return vz(); } float dxy() const { return dxy_; } float caloEta() const { return caloEta_; } diff --git a/DataFormats/L1TParticleFlow/interface/bit_encoding.h b/DataFormats/L1TParticleFlow/interface/bit_encoding.h index 9aa7446e75962..695594179d573 100644 --- a/DataFormats/L1TParticleFlow/interface/bit_encoding.h +++ b/DataFormats/L1TParticleFlow/interface/bit_encoding.h @@ -58,4 +58,34 @@ inline void l1pf_pattern_unpack(const ap_uint data[], T objs[N]) { } } +template +inline void l1pf_pattern_pack_slim(const T objs[N], ap_uint data[]) { +#ifdef __SYNTHESIS__ +#pragma HLS inline +#pragma HLS inline region recursive +#endif + assert(T::BITWIDTH_SLIM <= NB); + for (unsigned int i = 0; i < N; ++i) { +#ifdef __SYNTHESIS__ +#pragma HLS unroll +#endif + data[i + OFFS] = objs[i].pack_slim(); + } +} + +template +inline void l1pf_pattern_unpack_slim(const ap_uint data[], T objs[N]) { +#ifdef __SYNTHESIS__ +#pragma HLS inline +#pragma HLS inline region recursive +#endif + assert(T::BITWIDTH_SLIM <= NB); + for (unsigned int i = 0; i < N; ++i) { +#ifdef __SYNTHESIS__ +#pragma HLS unroll +#endif + objs[i] = T::unpack(data[i + OFFS]); + } +} + #endif diff --git a/DataFormats/L1TParticleFlow/interface/datatypes.h b/DataFormats/L1TParticleFlow/interface/datatypes.h index e579cd8493163..985e5ea6a4a9c 100644 --- a/DataFormats/L1TParticleFlow/interface/datatypes.h +++ b/DataFormats/L1TParticleFlow/interface/datatypes.h @@ -40,10 +40,10 @@ namespace l1ct { typedef ap_uint<13> tk2calo_dq_t; typedef ap_uint<4> egquality_t; typedef ap_uint<3> stub_t; - typedef ap_ufixed<12, 10, AP_TRN, AP_SAT> chi2_t; typedef ap_ufixed<10, 1, AP_TRN, AP_SAT> srrtot_t; - typedef ap_uint<8> meanz_t; // mean - SCALE_MEANZ = 320 + typedef ap_uint<8> meanz_t; // mean - MEANZ_OFFSET(= 320 cm) typedef ap_ufixed<10, 5, AP_TRN, AP_SAT> hoe_t; + typedef ap_uint<4> redChi2Bin_t; // FIXME: adjust range 10-11bits -> 1/4 - 1/2TeV is probably more than enough for all reasonable use cases typedef ap_ufixed<11, 9, AP_TRN, AP_SAT> iso_t; @@ -155,9 +155,9 @@ namespace l1ct { constexpr float Z0_LSB = 0.05; constexpr float DXY_LSB = 0.05; constexpr float PUPPIW_LSB = 1.0 / 256; - constexpr float MEANZ_SCALE = 320.; - constexpr float SRRTOT_LSB = pow(2,-9); - constexpr float HOE_LSB = pow(2,-5); + constexpr float MEANZ_OFFSET = 320.; + constexpr float SRRTOT_LSB = pow(2, -9); + constexpr float HOE_LSB = pow(2, -5); inline float floatPt(pt_t pt) { return pt.to_float(); } inline float floatPt(dpt_t pt) { return pt.to_float(); } @@ -174,9 +174,8 @@ namespace l1ct { inline float floatDxy(dxy_t dxy) { return dxy.to_float() * DXY_LSB; } inline float floatPuppiW(puppiWgt_t puppiw) { return puppiw.to_float() * PUPPIW_LSB; } inline float floatIso(iso_t iso) { return iso.to_float(); } - inline float floatChi2(chi2_t chi2) { return chi2.to_float(); } inline float floatSrrTot(srrtot_t srrtot) { return srrtot.to_float(); }; - inline float floatMeanZ(meanz_t meanz) { return meanz + MEANZ_SCALE; }; + inline float floatMeanZ(meanz_t meanz) { return meanz + MEANZ_OFFSET; }; inline float floatHoe(hoe_t hoe) { return hoe.to_float(); }; inline pt_t makePt(int pt) { return ap_ufixed<16, 14>(pt) >> 2; } @@ -208,10 +207,9 @@ namespace l1ct { inline iso_t makeIso(float iso) { return iso_t(0.25 * round(iso * 4)); } inline int makeDR2FromFloatDR(float dr) { return ceil(dr * dr / ETAPHI_LSB / ETAPHI_LSB); } - inline chi2_t makeChi2(float chi2) { return chi2_t(chi2); } - inline srrtot_t makeSrrTot(float var) { return srrtot_t(var); }; - inline meanz_t makeMeanZ(float var) { return round(var - MEANZ_SCALE); }; - inline hoe_t makeHoe(float var) { return hoe_t(var); }; + inline srrtot_t makeSrrTot(float var) { return srrtot_t(SRRTOT_LSB * round(var / SRRTOT_LSB)); }; + inline meanz_t makeMeanZ(float var) { return round(var - MEANZ_OFFSET); }; + inline hoe_t makeHoe(float var) { return hoe_t(HOE_LSB * round(var / HOE_LSB)); }; inline float maxAbsEta() { return ((1 << (eta_t::width - 1)) - 1) * ETAPHI_LSB; } inline float maxAbsPhi() { return ((1 << (phi_t::width - 1)) - 1) * ETAPHI_LSB; } diff --git a/DataFormats/L1TParticleFlow/interface/layer1_emulator.h b/DataFormats/L1TParticleFlow/interface/layer1_emulator.h index 152456059b939..944b2211e5f15 100644 --- a/DataFormats/L1TParticleFlow/interface/layer1_emulator.h +++ b/DataFormats/L1TParticleFlow/interface/layer1_emulator.h @@ -39,7 +39,7 @@ namespace l1ct { }; struct TkObjEmu : public TkObj { - // uint16_t hwChi2, hwStubs; + uint16_t hwChi2; float simPt, simCaloEta, simCaloPhi, simVtxEta, simVtxPhi, simZ0, simD0; const l1t::PFTrack *src = nullptr; bool read(std::fstream &from); @@ -47,8 +47,7 @@ namespace l1ct { void clear() { TkObj::clear(); src = nullptr; - // hwChi2 = 0; - // hwStubs = 0; + hwChi2 = 0; simPt = 0; simCaloEta = 0; simCaloPhi = 0; diff --git a/DataFormats/L1TParticleFlow/interface/layer1_objs.h b/DataFormats/L1TParticleFlow/interface/layer1_objs.h index a5aab55dcbf8f..28a008eba1207 100644 --- a/DataFormats/L1TParticleFlow/interface/layer1_objs.h +++ b/DataFormats/L1TParticleFlow/interface/layer1_objs.h @@ -49,8 +49,10 @@ namespace l1ct { bool hwIsEM() const { return hwEmID != 0; } - static const int BITWIDTH = pt_t::width + eta_t::width + phi_t::width + pt_t::width + emid_t::width + - srrtot_t::width + meanz_t::width + hoe_t::width; + static const int BITWIDTH_SLIM = pt_t::width + eta_t::width + phi_t::width + pt_t::width + emid_t::width; + + static const int BITWIDTH = BITWIDTH_SLIM + srrtot_t::width + meanz_t::width + hoe_t::width; + inline ap_uint pack() const { ap_uint ret; unsigned int start = 0; @@ -64,6 +66,7 @@ namespace l1ct { pack_into_bits(ret, start, hwHoe); return ret; } + inline static HadCaloObj unpack(const ap_uint &src) { HadCaloObj ret; unsigned int start = 0; @@ -77,6 +80,8 @@ namespace l1ct { unpack_from_bits(src, start, ret.hwHoe); return ret; } + + inline ap_uint pack_slim() const { return pack()(BITWIDTH_SLIM - 1, 0); } }; inline void clear(HadCaloObj &c) { c.clear(); } @@ -121,8 +126,10 @@ namespace l1ct { float floatMeanZ() const { return Scales::floatMeanZ(hwMeanZ); }; float floatHoe() const { return Scales::floatHoe(hwHoe); }; - static const int BITWIDTH = pt_t::width + eta_t::width + phi_t::width + pt_t::width + emid_t::width + - srrtot_t::width + meanz_t::width + hoe_t::width; + static const int BITWIDTH_SLIM = pt_t::width + eta_t::width + phi_t::width + pt_t::width + emid_t::width; + + static const int BITWIDTH = BITWIDTH_SLIM + srrtot_t::width + meanz_t::width + hoe_t::width; + inline ap_uint pack() const { ap_uint ret; unsigned int start = 0; @@ -150,6 +157,8 @@ namespace l1ct { return ret; } + + inline ap_uint pack_slim() const { return pack()(BITWIDTH_SLIM - 1, 0); } }; inline void clear(EmCaloObj &c) { c.clear(); } @@ -163,9 +172,11 @@ namespace l1ct { z0_t hwZ0; dxy_t hwDxy; tkquality_t hwQuality; - // FIXME: these variables are actually not in the track word or not in this format... stub_t hwStubs; - chi2_t hwChi2; + redChi2Bin_t hwRedChi2RZ; // 4 bits + redChi2Bin_t hwRedChi2RPhi; // 4 bits + //FIXME: is this actually filled? 3 bits would be enough + redChi2Bin_t hwRedChi2Bend; // 4 bits enum TkQuality { PFLOOSE = 1, PFTIGHT = 2 }; bool isPFLoose() const { return hwQuality[0]; } @@ -176,7 +187,8 @@ namespace l1ct { inline bool operator==(const TkObj &other) const { return hwPt == other.hwPt && hwEta == other.hwEta && hwPhi == other.hwPhi && hwDEta == other.hwDEta && hwDPhi == other.hwDPhi && hwZ0 == other.hwZ0 && hwDxy == other.hwDxy && hwCharge == other.hwCharge && - hwQuality == other.hwQuality && hwChi2 == other.hwChi2 && hwStubs == other.hwStubs; + hwQuality == other.hwQuality && hwStubs == other.hwStubs && hwRedChi2RZ == other.hwRedChi2RZ && + hwRedChi2RPhi == other.hwRedChi2RPhi && hwRedChi2Bend == other.hwRedChi2Bend; } inline bool operator>(const TkObj &other) const { return hwPt > other.hwPt; } @@ -193,7 +205,9 @@ namespace l1ct { hwCharge = false; hwQuality = 0; hwStubs = 0; - hwChi2 = 0; + hwRedChi2RZ = 0; + hwRedChi2RPhi = 0; + hwRedChi2Bend = 0; } int intPt() const { return Scales::intPt(hwPt); } @@ -211,10 +225,13 @@ namespace l1ct { float floatVtxPhi() const { return Scales::floatPhi(hwVtxPhi()); } float floatZ0() const { return Scales::floatZ0(hwZ0); } float floatDxy() const { return Scales::floatDxy(hwDxy); } - float floatChi2() const { return Scales::floatChi2(hwChi2); } - static const int BITWIDTH = pt_t::width + eta_t::width + phi_t::width + tkdeta_t::width + tkdphi_t::width + 1 + - z0_t::width + dxy_t::width + tkquality_t::width + stub_t::width + chi2_t::width; + static const int BITWIDTH_SLIM = pt_t::width + eta_t::width + phi_t::width + tkdeta_t::width + tkdphi_t::width + 1 + + z0_t::width + dxy_t::width + tkquality_t::width; + + static const int BITWIDTH = + BITWIDTH_SLIM + stub_t::width + redChi2Bin_t::width + redChi2Bin_t::width + redChi2Bin_t::width; + inline ap_uint pack() const { ap_uint ret; unsigned int start = 0; @@ -228,7 +245,9 @@ namespace l1ct { pack_into_bits(ret, start, hwDxy); pack_into_bits(ret, start, hwQuality); pack_into_bits(ret, start, hwStubs); - pack_into_bits(ret, start, hwChi2); + pack_into_bits(ret, start, hwRedChi2RZ); + pack_into_bits(ret, start, hwRedChi2RPhi); + pack_into_bits(ret, start, hwRedChi2Bend); return ret; } inline static TkObj unpack(const ap_uint &src) { @@ -244,9 +263,12 @@ namespace l1ct { unpack_from_bits(src, start, ret.hwDxy); unpack_from_bits(src, start, ret.hwQuality); unpack_from_bits(src, start, ret.hwStubs); - unpack_from_bits(src, start, ret.hwChi2); + unpack_from_bits(src, start, ret.hwRedChi2RZ); + unpack_from_bits(src, start, ret.hwRedChi2RPhi); + unpack_from_bits(src, start, ret.hwRedChi2Bend); return ret; } + inline ap_uint pack_slim() const { return pack()(BITWIDTH_SLIM - 1, 0); } }; inline void clear(TkObj &c) { c.clear(); } diff --git a/DataFormats/L1TParticleFlow/src/layer1_emulator.cpp b/DataFormats/L1TParticleFlow/src/layer1_emulator.cpp index f0b56b1c25022..3f1284413bbf7 100644 --- a/DataFormats/L1TParticleFlow/src/layer1_emulator.cpp +++ b/DataFormats/L1TParticleFlow/src/layer1_emulator.cpp @@ -34,14 +34,14 @@ bool l1ct::EmCaloObjEmu::write(std::fstream& to) const { return writeObj(from, *this) && readVar(from, hwChi2) && readVar(from, hwStubs) && readVar(from, simPt) && - readVar(from, simCaloEta) && readVar(from, simCaloPhi) && readVar(from, simVtxEta) && - readVar(from, simVtxPhi) && readVar(from, simZ0) && readVar(from, simD0); + return readObj(from, *this) && readVar(from, hwChi2) && readVar(from, simPt) && readVar(from, simCaloEta) && + readVar(from, simCaloPhi) && readVar(from, simVtxEta) && readVar(from, simVtxPhi) && readVar(from, simZ0) && + readVar(from, simD0); } bool l1ct::TkObjEmu::write(std::fstream& to) const { - return writeObj(*this, to) && writeVar(hwChi2, to) && writeVar(hwStubs, to) && writeVar(simPt, to) && - writeVar(simCaloEta, to) && writeVar(simCaloPhi, to) && writeVar(simVtxEta, to) && writeVar(simVtxPhi, to) && - writeVar(simZ0, to) && writeVar(simD0, to); + return writeObj(*this, to) && writeVar(hwChi2, to) && writeVar(simPt, to) && writeVar(simCaloEta, to) && + writeVar(simCaloPhi, to) && writeVar(simVtxEta, to) && writeVar(simVtxPhi, to) && writeVar(simZ0, to) && + writeVar(simD0, to); } bool l1ct::MuObjEmu::read(std::fstream& from) { diff --git a/L1Trigger/L1THGCalUtilities/test/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc b/L1Trigger/L1THGCalUtilities/test/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc index 6a495b3e5e96e..34c72aa4df886 100644 --- a/L1Trigger/L1THGCalUtilities/test/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc +++ b/L1Trigger/L1THGCalUtilities/test/ntuples/HGCalTriggerNtupleHGCMulticlusters.cc @@ -181,12 +181,12 @@ void HGCalTriggerNtupleHGCMulticlusters::fill(const edm::Event& e, const HGCalTr cl3d_layer_pt_.emplace_back(layer_pt); } - // // Retrieve indices of trigger cells inside cluster - // cl3d_clusters_id_.emplace_back(cl3d_itr->constituents().size()); - // std::transform(cl3d_itr->constituents_begin(), - // cl3d_itr->constituents_end(), - // cl3d_clusters_id_.back().begin(), - // [](const std::pair>& id_cl) { return id_cl.second->detId(); }); + // Retrieve indices of trigger cells inside cluster + cl3d_clusters_id_.emplace_back(cl3d_itr->constituents().size()); + std::transform(cl3d_itr->constituents_begin(), + cl3d_itr->constituents_end(), + cl3d_clusters_id_.back().begin(), + [](const std::pair>& id_cl) { return id_cl.second->detId(); }); } } diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h b/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h index 58a8492bd4def..ab91b6ff28fed 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h @@ -1,61 +1,60 @@ #ifndef CC_INVERSION_H__ #define CC_INVERSION_H__ -constexpr int ceillog2(int x){ - return (x <= 2) ? 1 : 1 + ceillog2((x+1) / 2); -} +constexpr int ceillog2(int x) { return (x <= 2) ? 1 : 1 + ceillog2((x + 1) / 2); } -template -inline float real_val_from_idx(unsigned i){ - // Treat the index as the top N bits - static constexpr int NB = ceillog2(N); // number of address bits for table - data_T x(0); - // The MSB of 1 is implicit in the table - x[x.width-1] = 1; - // So we can use the next NB bits for real data - x(x.width-2, x.width-NB-1) = i; - return (float) x; +template +inline float real_val_from_idx(unsigned i) { + // Treat the index as the top N bits + static constexpr int NB = ceillog2(N); // number of address bits for table + data_T x(0); + // The MSB of 1 is implicit in the table + x[x.width - 1] = 1; + // So we can use the next NB bits for real data + x(x.width - 2, x.width - NB - 1) = i; + return (float)x; } -template -inline unsigned idx_from_real_val(data_T x){ - // Slice the top N bits to get an index into the table - static constexpr int NB = ceillog2(N); // number of address bits for table - // Slice the top-1 NB bits of the value - // the MSB of '1' is implicit, so only slice below that - ap_uint y = x(x.width-2, x.width-NB-1); - return (unsigned) y(NB-1, 0); +template +inline unsigned idx_from_real_val(data_T x) { + // Slice the top N bits to get an index into the table + static constexpr int NB = ceillog2(N); // number of address bits for table + // Slice the top-1 NB bits of the value + // the MSB of '1' is implicit, so only slice below that + ap_uint y = x(x.width - 2, x.width - NB - 1); + return (unsigned)y(NB - 1, 0); } -template -void init_invert_table(table_T table_out[N]){ - // The template data_T is the data type used to address the table - for(unsigned i = 0; i < N; i++){ - float x = real_val_from_idx(i); - table_T inv_x = 1 / x; - table_out[i] = inv_x; - } +template +void init_invert_table(table_T table_out[N]) { + // The template data_T is the data type used to address the table + for (unsigned i = 0; i < N; i++) { + float x = real_val_from_idx(i); + table_T inv_x = 1 / x; + table_out[i] = inv_x; + } } -template -table_t invert_with_shift(in_t in){ - table_t inv_table[N]; - init_invert_table(inv_table); +template +table_t invert_with_shift(in_t in) { + table_t inv_table[N]; + init_invert_table(inv_table); - // find the first '1' in the denominator - int msb = 0; - for(int b = 0; b < in.width; b++){ - #pragma HLS unroll - if(in[b]) msb = b; - } - // shift up the denominator such that the left-most bit (msb) is '1' - in_t in_shifted = in << (in.width-msb-1); - // lookup the inverse of the shifted input - int idx = idx_from_real_val(in_shifted); - table_t inv_in = inv_table[idx]; - // shift the output back - table_t out = inv_in << (in.width-msb-1); - return out; + // find the first '1' in the denominator + int msb = 0; + for (int b = 0; b < in.width; b++) { +// #pragma HLS unroll + if (in[b]) + msb = b; + } + // shift up the denominator such that the left-most bit (msb) is '1' + in_t in_shifted = in << (in.width - msb - 1); + // lookup the inverse of the shifted input + int idx = idx_from_real_val(in_shifted); + table_t inv_in = inv_table[idx]; + // shift the output back + table_t out = inv_in << (in.width - msb - 1); + return out; } #endif \ No newline at end of file diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/conifer.h b/L1Trigger/Phase2L1ParticleFlow/interface/conifer.h index 5b53956a4d032..f2c6efedc4676 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/conifer.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/conifer.h @@ -1,8 +1,12 @@ #ifndef L1TRIGGER_PHASE2L1PARTICLEFLOW_CONNIFER_H #define L1TRIGGER_PHASE2L1PARTICLEFLOW_CONNIFER_H -//#include "FWCore/Utilities/interface/Exception.h" #include "nlohmann/json.hpp" #include +#ifdef CMSSW_GIT_HASH +#include "FWCore/Utilities/interface/Exception.h" +#else +#include +#endif namespace conifer { @@ -115,18 +119,25 @@ namespace conifer { std::vector decision_function(std::vector x) const { /* Do the prediction */ - //if (x.size() != n_features) { - // throw cms::Exception("RuntimeError") - // << "Conifer : Size of feature vector mismatches expected n_features" << std::endl; - //} +#ifdef CMSSW_GIT_HASH + if (x.size() != n_features) { + throw cms::Exception("RuntimeError") + << "Conifer : Size of feature vector mismatches expected n_features" << std::endl; + } +#else + if (x.size() != n_features) { + throw std::runtime_error("Conifer : Size of feature vector mismatches expected n_features"); + } +#endif std::vector values; std::vector> values_trees; values_trees.resize(n_classes); values.resize(n_classes, U(0)); for (unsigned int i = 0; i < n_classes; i++) { - std::transform(trees.begin(), trees.end(), std::back_inserter(values_trees.at(i)), [&i, &x](std::vector> tree_v) { - return tree_v.at(i).decision_function(x); - }); + std::transform(trees.begin(), + trees.end(), + std::back_inserter(values_trees.at(i)), + [&i, &x](std::vector> tree_v) { return tree_v.at(i).decision_function(x); }); if (useAddTree) { values.at(i) = init_predict_.at(i); values.at(i) += reduce>(values_trees.at(i), add); diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h index 7eb6cf5e449f8..521722b9c5a19 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h @@ -33,7 +33,7 @@ namespace l1ct { std::vector dPhiValues; float trkQualityPtMin; // GeV bool doCompositeTkEle; - unsigned int nCOMPCAND_PER_CLUSTER; + unsigned int nCompCandPerCluster; bool writeEgSta; struct IsoParameters { @@ -58,43 +58,13 @@ namespace l1ct { EGIsoEleObjEmu::IsoType hwIsoTypeTkEle; EGIsoObjEmu::IsoType hwIsoTypeTkEm; - //bool doCompositeTkEle; struct CompIDParameters { CompIDParameters(const edm::ParameterSet &); - CompIDParameters(double hoeMin, double hoeMax, double tkptMin, double tkptMax, double srrtotMin, double srrtotMax, double detaMin, double detaMax, double dptMin, double dptMax, double meanzMin, double meanzMax, double dphiMin, double dphiMax, double tkchi2Min, double tkchi2Max, double tkz0Min, double tkz0Max, double tknstubsMin, double tknstubsMax, double BDTcut_wp97p5, double BDTcut_wp95p0) - : hoeMin(hoeMin), hoeMax(hoeMax), - tkptMin(tkptMin),tkptMax(tkptMax), - srrtotMin(srrtotMin),srrtotMax(srrtotMax), - detaMin(detaMin),detaMax(detaMax), - dptMin(dptMin),dptMax(dptMax), - meanzMin(meanzMin),meanzMax(meanzMax), - dphiMin(dphiMin),dphiMax(dphiMax), - tkchi2Min(tkchi2Min),tkchi2Max(tkchi2Max), - tkz0Min(tkz0Min),tkz0Max(tkz0Max), - tknstubsMin(tknstubsMin),tknstubsMax(tknstubsMax), - BDTcut_wp97p5(BDTcut_wp97p5),BDTcut_wp95p0(BDTcut_wp95p0){} - double hoeMin; - double hoeMax; - double tkptMin; - double tkptMax; - double srrtotMin; - double srrtotMax; - double detaMin; - double detaMax; - double dptMin; - double dptMax; - double meanzMin; - double meanzMax; - double dphiMin; - double dphiMax; - double tkchi2Min; - double tkchi2Max; - double tkz0Min; - double tkz0Max; - double tknstubsMin; - double tknstubsMax; - double BDTcut_wp97p5; - double BDTcut_wp95p0; + CompIDParameters(double bdtScore_loose_wp, double bdtScore_tight_wp, const std::string &model) + : bdtScore_loose_wp(bdtScore_loose_wp), bdtScore_tight_wp(bdtScore_tight_wp), conifer_model(model) {} + const double bdtScore_loose_wp; // XGBOOST score + const double bdtScore_tight_wp; // XGBOOST score + const std::string conifer_model; }; CompIDParameters compIDparams; @@ -130,7 +100,7 @@ namespace l1ct { EGIsoEleObjEmu::IsoType hwIsoTypeTkEle = EGIsoEleObjEmu::IsoType::TkIso, EGIsoObjEmu::IsoType hwIsoTypeTkEm = EGIsoObjEmu::IsoType::TkIsoPV, // FIXME: maybe we round these? - const CompIDParameters &myCompIDparams = {-1.0, 1566.547607421875, 1.9501149654388428, 11102.0048828125, 0.0, 0.01274710614234209, -0.24224889278411865, 0.23079538345336914, 0.010325592942535877, 184.92538452148438, 325.0653991699219, 499.6089782714844, -6.281332015991211, 6.280326843261719, 0.024048099294304848, 1258.37158203125, -14.94140625, 14.94140625, 4.0, 6.0, 0.7927004, 0.9826955}, + const CompIDParameters &compIDparams = {0.7927004, 0.9826955, "compositeID.json"}, int debug = 0) : nTRACK(nTrack), @@ -153,7 +123,7 @@ namespace l1ct { dPhiValues(dPhiValues), trkQualityPtMin(trkQualityPtMin), doCompositeTkEle(doCompositeTkEle), - nCOMPCAND_PER_CLUSTER(nCompCandPerCluster), + nCompCandPerCluster(nCompCandPerCluster), writeEgSta(writeEgSta), tkIsoParams_tkEle(tkIsoParams_tkEle), tkIsoParams_tkEm(tkIsoParams_tkEm), @@ -163,16 +133,14 @@ namespace l1ct { doPfIso(doPfIso), hwIsoTypeTkEle(hwIsoTypeTkEle), hwIsoTypeTkEm(hwIsoTypeTkEm), - compIDparams(myCompIDparams), + compIDparams(compIDparams), debug(debug) {} }; class PFTkEGAlgoEmulator { public: - PFTkEGAlgoEmulator(const PFTkEGAlgoEmuConfig &config); - virtual ~PFTkEGAlgoEmulator() {} void toFirmware(const PFInputRegion &in, PFRegion ®ion, EmCaloObj calo[/*nCALO*/], TkObj track[/*nTRACK*/]) const; @@ -190,26 +158,26 @@ namespace l1ct { bool writeEgSta() const { return cfg.writeEgSta; } - private: - + typedef ap_fixed<21, 12, AP_RND_CONV, AP_SAT> bdt_feature_t; + private: void link_emCalo2emCalo(const std::vector &emcalo, std::vector &emCalo2emCalo) const; void link_emCalo2tk_elliptic(const PFRegionEmu &r, - const std::vector &emcalo, - const std::vector &track, - std::vector &emCalo2tk) const; + const std::vector &emcalo, + const std::vector &track, + std::vector &emCalo2tk) const; void link_emCalo2tk_composite(const PFRegionEmu &r, - const std::vector &emcalo, - const std::vector &track, - std::vector &emCalo2tk, - std::vector &emCaloTkBdtScore) const; + const std::vector &emcalo, + const std::vector &track, + std::vector &emCalo2tk, + std::vector &emCaloTkBdtScore) const; struct CompositeCandidate { unsigned int cluster_idx; unsigned int track_idx; - double dpt; // For sorting + double dpt; // For sorting }; float compute_composite_score(CompositeCandidate &cand, @@ -389,7 +357,7 @@ namespace l1ct { z0_t z0) const; PFTkEGAlgoEmuConfig cfg; - conifer::BDT,ap_fixed<12,3,AP_RND_CONV,AP_SAT>,0> * composite_bdt_; + conifer::BDT, false> *composite_bdt_; int debug_; }; } // namespace l1ct diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/tkinput_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/tkinput_ref.h index 585156036ca9a..2545a6dabc243 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/tkinput_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/tkinput_ref.h @@ -22,14 +22,24 @@ namespace l1ct { }; TrackInputEmulator(const edm::ParameterSet &iConfig); - TrackInputEmulator(Region = Region::Endcap, Encoding encoding = Encoding::Stepping, bool bitwise = true); + TrackInputEmulator(Region = Region::Endcap, + Encoding encoding = Encoding::Stepping, + bool bitwise = true, + bool slim = true); std::pair decodeTrack(ap_uint<96> tkword, const l1ct::PFRegionEmu §or) const { - return decodeTrack(tkword, sector, bitwise_); + return decodeTrack(tkword, sector, bitwise_, slim_); } std::pair decodeTrack(ap_uint<96> tkword, const l1ct::PFRegionEmu §or, - bool bitwise) const; + bool bitwise) const { + return decodeTrack(tkword, sector, bitwise, slim_); + } + + std::pair decodeTrack(ap_uint<96> tkword, + const l1ct::PFRegionEmu §or, + bool bitwise, + bool slim) const; //== Unpackers == static bool valid(const ap_uint<96> &tkword) { return tkword[95]; } @@ -151,6 +161,15 @@ namespace l1ct { const std::vector &tanlLUT() const { return tanlLUT_; } const std::vector &ptLUT() const { return ptLUT_; } + unsigned int countSetBits(unsigned int n) const { + unsigned int count = 0; + while (n) { + n &= (n - 1); + count++; + } + return count; + } + protected: // utilities template @@ -175,6 +194,9 @@ namespace l1ct { /// Whether to run the bitwise accurate or floating point conversions bool bitwise_; + /// Whether to unpack and populate also nstubs and various chi2 variables (needed for CompibedID in the endcap) + bool slim_; + /// Main constants float rInvToPt_, phiScale_, z0Scale_; diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc index 3031ffd06460f..15d0a72e83b50 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc @@ -176,9 +176,7 @@ L1TCorrelatorLayer1Producer::L1TCorrelatorLayer1Producer(const edm::ParameterSet #if 0 // LATER produces("TKVtx"); #endif -#if 0 // LATER produces>("DecodedTK"); -#endif for (const auto &tag : iConfig.getParameter>("emClusters")) { emCands_.push_back(consumes(tag)); @@ -368,10 +366,7 @@ void L1TCorrelatorLayer1Producer::produce(edm::Event &iEvent, const edm::EventSe iEvent.put(fetchEmCalo(), "EmCalo"); iEvent.put(fetchHadCalo(), "Calo"); iEvent.put(fetchTracks(), "TK"); - - #if 0 - iEvent.put(fetchDecodedTracks(), "DecodedTK"); - #endif + iEvent.put(fetchDecodedTracks(), "DecodedTK"); // Then do the vertexing, and save it out std::vector z0s; @@ -643,8 +638,7 @@ void L1TCorrelatorLayer1Producer::addDecodedTrack(l1ct::DetectorSector 0; } // CMSSW-only extra info - tkAndSel.first.hwChi2 = l1ct::Scales::makeChi2(t.chi2()); - tkAndSel.first.hwStubs = t.nStubs(); + tkAndSel.first.hwChi2 = round(t.chi2() * 10); tkAndSel.first.simPt = t.pt(); tkAndSel.first.simCaloEta = t.caloEta(); tkAndSel.first.simCaloPhi = t.caloPhi(); @@ -653,7 +647,7 @@ void L1TCorrelatorLayer1Producer::addDecodedTrack(l1ct::DetectorSector> &sec, const l1t::PFCluster &c) { ap_uint<256> cwrd = 0; - ap_uint<14> w_pt = round(c.pt() / 0.25); + ap_ufixed<14, 12, AP_RND_CONV, AP_SAT> w_pt = c.pt(); ap_uint<14> w_empt = round(c.emEt() / 0.25); constexpr float ETAPHI_LSB = M_PI / 720; ap_int<9> w_eta = round(sec.region.localEta(c.eta()) / ETAPHI_LSB); ap_int<9> w_phi = round(sec.region.localPhi(c.phi()) / ETAPHI_LSB); ap_uint<10> w_qual = c.hwQual(); - ap_uint<13> w_srrtot = round(c.sigmaRR()/l1ct::Scales::SRRTOT_LSB); + ap_uint<13> w_srrtot = round(c.sigmaRR() / l1ct::Scales::SRRTOT_LSB); ap_uint<12> w_meanz = round(c.absZBarycenter()); - ap_uint<12> w_hoe = round(c.hOverE()/l1ct::Scales::HOE_LSB); - - cwrd(13, 0) = w_pt; + // FIXME: the calibration can actually make hoe become negative....we add a small protection for now + // We use ap_ufixed to handle saturation and rounding + ap_ufixed<10, 5, AP_RND_CONV, AP_SAT> w_hoe = c.hOverE(); + + cwrd(13, 0) = w_pt.range(); cwrd(27, 14) = w_empt; cwrd(72, 64) = w_eta; cwrd(81, 73) = w_phi; cwrd(115, 106) = w_qual; - // FIXME: we add the variables use by composite-ID. The definitin will have to be reviewd once the + // FIXME: we add the variables use by composite-ID. The definitin will have to be reviewd once the // hgc format is better defined. For now we use // hwMeanZ = word 1 bits 30-19 // hwSrrTot = word 3 bits 21 - 9 - // hoe = word 1 bits 63-52 (currently spare) + // hoe = word 1 bits 63-52 (currently spare in the interface) cwrd(213, 201) = w_srrtot; cwrd(94, 83) = w_meanz; - // FIXME: we use a spare space in the word for hoe which is not in the current interface - cwrd(127, 116) = w_hoe; + cwrd(127, 116) = w_hoe.range(); sec.obj.push_back(cwrd); } @@ -871,15 +866,16 @@ std::unique_ptr L1TCorrelatorLayer1Producer::fetchTr std::unique_ptr> L1TCorrelatorLayer1Producer::fetchDecodedTracks() const { auto ret = std::make_unique>(); - for (const auto r : event_.decoded.track) { + for (const auto &r : event_.decoded.track) { const auto ® = r.region; for (const auto &p : r.obj) { if (p.hwPt == 0 || !reg.isFiducial(p)) continue; - reco::Particle::PolarLorentzVector p4(p.floatPt(), reg.floatGlbEta(p.hwVtxEta()), reg.floatGlbPhi(p.hwVtxPhi()), 0); + reco::Particle::PolarLorentzVector p4( + p.floatPt(), reg.floatGlbEta(p.hwVtxEta()), reg.floatGlbPhi(p.hwVtxPhi()), 0); + + reco::Particle::Point vtx(0, 0, p.floatZ0()); - reco::Particle::Point vtx(0,0,p.floatZ0()); - ret->emplace_back(l1t::PFTrack(p.intCharge(), reco::Particle::LorentzVector(p4), vtx, @@ -899,7 +895,6 @@ std::unique_ptr> L1TCorrelatorLayer1Producer::fetchDec return ret; } - std::unique_ptr L1TCorrelatorLayer1Producer::fetchPF() const { auto ret = std::make_unique(); for (unsigned int ir = 0, nr = event_.pfinputs.size(); ir < nr; ++ir) { @@ -922,7 +917,7 @@ std::unique_ptr L1TCorrelatorLayer1Producer::fetchPF ret->back().setHwTkQuality(p.hwTkQuality); ret->back().setCaloEta(reg.floatGlbEtaOf(p)); ret->back().setCaloPhi(reg.floatGlbPhiOf(p)); - + setRefs_(ret->back(), p); } for (const auto &p : event_.out[ir].pfneutral) { diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/PFClusterProducerFromHGC3DClusters.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/PFClusterProducerFromHGC3DClusters.cc index 8594fc19f3c16..0a76d751f90c0 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/PFClusterProducerFromHGC3DClusters.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/PFClusterProducerFromHGC3DClusters.cc @@ -129,7 +129,8 @@ void l1tpf::PFClusterProducerFromHGC3DClusters::produce(edm::Event &iEvent, cons //float em_old = cluster.emEt(); float em_new = it->iPt(l1t::HGCalMulticluster::EnergyInterpretation::EM); float pt_new = had_old + em_new; - float hoe_new = em_new > 0 ? (had_old / em_new) : -1; + // FIXME: -1 can be a problem for later stages of the processing. For now we set it to something which saturates the hoe variable + float hoe_new = em_new > 0 ? (had_old / em_new) : 999; cluster = l1t::PFCluster(pt_new, it->eta(), it->phi(), hoe_new, /*isEM=*/isEM); //printf("Scenario %d: pt %7.2f eta %+5.3f em %7.2f, EMI %7.2f, h/e % 8.3f --> pt %7.2f, em %7.2f, h/e % 8.3f\n", // 2, pt, it->eta(), em_old, em_new, hoe, cluster.pt(), cluster.emEt(), cluster.hOverE()); diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py index e270bb99bca73..96693cc0f86ae 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py @@ -1,5 +1,8 @@ import FWCore.ParameterSet.Config as cms +import math +def BDTScore_tmva2xgboost(score): + return 1. / (1. + math.sqrt((1. - score) / (1. + score))) tkEgAlgoParameters = cms.PSet( nTRACK=cms.uint32(50), # very large numbers for first test @@ -52,30 +55,11 @@ hwIsoTypeTkEle=cms.uint32(0), hwIsoTypeTkEm=cms.uint32(2), doCompositeTkEle=cms.bool(False), - nCOMPCAND_PER_CLUSTER=cms.uint32(3), - compositeParametersTkEle=cms.PSet( # Parameters used to normalize input features - hoeMin=cms.double(-1.0), - hoeMax=cms.double(1566.547607421875), - tkptMin=cms.double(1.9501149654388428), - tkptMax=cms.double(11102.0048828125), - srrtotMin=cms.double(0.0), - srrtotMax=cms.double(0.01274710614234209), - detaMin=cms.double(-0.24224889278411865), - detaMax=cms.double(0.23079538345336914), - dptMin=cms.double(0.010325592942535877), - dptMax=cms.double(184.92538452148438), - meanzMin=cms.double(325.0653991699219), - meanzMax=cms.double(499.6089782714844), - dphiMin=cms.double(-6.281332015991211), - dphiMax=cms.double(6.280326843261719), - tkchi2Min=cms.double(0.024048099294304848), - tkchi2Max=cms.double(1258.37158203125), - tkz0Min=cms.double(-14.94140625), - tkz0Max=cms.double(14.94140625), - tknstubsMin=cms.double(4.0), - tknstubsMax=cms.double(6.0), - BDTcut_wp97p5=cms.double(0.7927004), - BDTcut_wp95p0=cms.double(0.9826955), + nCompCandPerCluster=cms.uint32(3), + compositeParametersTkEle=cms.PSet( + bdt_loose_wp=cms.double(0.), + bdt_tight_wp=cms.double(BDTScore_tmva2xgboost(0.9826955)), + conifer_model=cms.string("L1Trigger/Phase2L1ParticleFlow/data/compositeID.json") ), ) diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgo_cfi.py b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgo_cfi.py deleted file mode 100644 index 500121131a360..0000000000000 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgo_cfi.py +++ /dev/null @@ -1,47 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -tkEgConfig = cms.PSet( - debug=cms.untracked.int32(0), - doBremRecovery=cms.bool(True), - doTkIsolation=cms.bool(True), - filterHwQuality=cms.bool(True), - caloHwQual=cms.int32(4), - dEtaMaxBrem=cms.double(0.02), - dPhiMaxBrem=cms.double(0.1), - absEtaBoundaries=cms.vdouble(0.0, 0.9, 1.5), - dEtaValues=cms.vdouble(0.025, 0.015, 0.01), # last was 0.0075 in TDR - dPhiValues=cms.vdouble(0.07, 0.07, 0.07), - caloEtMin=cms.double(0.0), - trkQualityPtMin=cms.double(10.0), - trkQualityChi2=cms.double(1e10), - writeEgSta=cms.bool(True), - tkIsoParametersTkEm=cms.PSet( - tkQualityPtMin=cms.double(2.), - dZ=cms.double(0.6), - dRMin=cms.double(0.07), - dRMax=cms.double(0.30), - tkQualityChi2Max=cms.double(100), - ), - tkIsoParametersTkEle=cms.PSet( - tkQualityPtMin=cms.double(2.), - dZ=cms.double(0.6), - dRMin=cms.double(0.03), - dRMax=cms.double(0.20), - tkQualityChi2Max=cms.double(1e10), - ), - pfIsoParametersTkEm=cms.PSet( - tkQualityPtMin=cms.double(1.), - dZ=cms.double(0.6), - dRMin=cms.double(0.07), - dRMax=cms.double(0.30), - tkQualityChi2Max=cms.double(100), - ), - pfIsoParametersTkEle=cms.PSet( - tkQualityPtMin=cms.double(1.), - dZ=cms.double(0.6), - dRMin=cms.double(0.03), - dRMax=cms.double(0.20), - tkQualityChi2Max=cms.double(1e10), - ) - -) diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py index abcf1200d706f..c1e1acd8cff18 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py @@ -29,6 +29,7 @@ region = cms.string("barrel"), trackWordEncoding = cms.string("biased"), bitwiseAccurate = cms.bool(True), + slimDataFormat = cms.bool(True), ptLUTBits = cms.uint32(11), etaLUTBits = cms.uint32(10), etaPreOffs = cms.int32(0), @@ -169,6 +170,7 @@ region = cms.string("endcap"), trackWordEncoding = cms.string("biased"), bitwiseAccurate = cms.bool(True), + slimDataFormat = cms.bool(False), ptLUTBits = cms.uint32(11), etaLUTBits = cms.uint32(11), etaPreOffs = cms.int32(0), diff --git a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp index 5786ac32bebe0..67d43529efcdd 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp @@ -9,10 +9,6 @@ #include #include -#include "DataFormats/L1TParticleFlow/interface/PFTrack.h" -#include "DataFormats/L1TParticleFlow/interface/PFCluster.h" - - using namespace l1ct; #ifdef CMSSW_GIT_HASH @@ -36,7 +32,7 @@ l1ct::PFTkEGAlgoEmuConfig::PFTkEGAlgoEmuConfig(const edm::ParameterSet &pset) dPhiValues(pset.getParameter>("dPhiValues")), trkQualityPtMin(pset.getParameter("trkQualityPtMin")), doCompositeTkEle(pset.getParameter("doCompositeTkEle")), - nCOMPCAND_PER_CLUSTER(pset.getParameter("nCOMPCAND_PER_CLUSTER")), + nCompCandPerCluster(pset.getParameter("nCompCandPerCluster")), writeEgSta(pset.getParameter("writeEGSta")), tkIsoParams_tkEle(pset.getParameter("tkIsoParametersTkEle")), tkIsoParams_tkEm(pset.getParameter("tkIsoParametersTkEm")), @@ -56,44 +52,21 @@ l1ct::PFTkEGAlgoEmuConfig::IsoParameters::IsoParameters(const edm::ParameterSet pset.getParameter("dRMax")) {} l1ct::PFTkEGAlgoEmuConfig::CompIDParameters::CompIDParameters(const edm::ParameterSet &pset) - : CompIDParameters(pset.getParameter("hoeMin"), - pset.getParameter("hoeMax"), - pset.getParameter("tkptMin"), - pset.getParameter("tkptMax"), - pset.getParameter("srrtotMin"), - pset.getParameter("srrtotMax"), - pset.getParameter("detaMin"), - pset.getParameter("detaMax"), - pset.getParameter("dptMin"), - pset.getParameter("dptMax"), - pset.getParameter("meanzMin"), - pset.getParameter("meanzMax"), - pset.getParameter("dphiMin"), - pset.getParameter("dphiMax"), - pset.getParameter("tkchi2Min"), - pset.getParameter("tkchi2Max"), - pset.getParameter("tkz0Min"), - pset.getParameter("tkz0Max"), - pset.getParameter("tknstubsMin"), - pset.getParameter("tknstubsMax"), - pset.getParameter("BDTcut_wp97p5"), - pset.getParameter("BDTcut_wp95p0")) {} + : CompIDParameters(pset.getParameter("bdt_loose_wp"), + pset.getParameter("bdt_tight_wp"), + pset.getParameter("conifer_model")) {} #endif -PFTkEGAlgoEmulator::PFTkEGAlgoEmulator(const PFTkEGAlgoEmuConfig &config) : cfg(config), -composite_bdt_(nullptr), -debug_(cfg.debug) { - if(cfg.doCompositeTkEle) { - //FIXME: make the name of the file configurable +PFTkEGAlgoEmulator::PFTkEGAlgoEmulator(const PFTkEGAlgoEmuConfig &config) + : cfg(config), composite_bdt_(nullptr), debug_(cfg.debug) { + if (cfg.doCompositeTkEle) { #ifdef CMSSW_GIT_HASH - auto resolvedFileName = edm::FileInPath("L1Trigger/Phase2L1ParticleFlow/data/compositeID.json").fullPath(); + auto resolvedFileName = edm::FileInPath(cfg.compIDparams.conifer_model).fullPath(); #else - auto resolvedFileName = "compositeID.json"; + auto resolvedFileName = cfg.compIDparams.conifer_model; #endif - std::cout<,ap_fixed<12,3,AP_RND_CONV,AP_SAT>,0> (resolvedFileName); - std::cout<<"declared bdt"<, false>(resolvedFileName); } } @@ -157,7 +130,6 @@ void PFTkEGAlgoEmulator::link_emCalo2emCalo(const std::vector &emc } } - void PFTkEGAlgoEmulator::link_emCalo2tk_elliptic(const PFRegionEmu &r, const std::vector &emcalo, const std::vector &track, @@ -195,70 +167,61 @@ void PFTkEGAlgoEmulator::link_emCalo2tk_elliptic(const PFRegionEmu &r, } } - void PFTkEGAlgoEmulator::link_emCalo2tk_composite(const PFRegionEmu &r, - const std::vector &emcalo, - const std::vector &track, - std::vector &emCalo2tk, - std::vector &emCaloTkBdtScore) const { + const std::vector &emcalo, + const std::vector &track, + std::vector &emCalo2tk, + std::vector &emCaloTkBdtScore) const { unsigned int nTrackMax = std::min(track.size(), cfg.nTRACK_EGIN); - std::cout<<"doing loose dR matching"< candidates; for (unsigned int itk = 0; itk < nTrackMax; ++itk) { - std::cout<<"track "< bool - { return a.dpt < b.dpt; }); - unsigned int nCandPerCluster = std::min(candidates.size(), cfg.nCOMPCAND_PER_CLUSTER); - std::cout << "# composite candidates: " << nCandPerCluster << std::endl; - if(nCandPerCluster == 0) continue; - - float bdtWP_MVA = cfg.compIDparams.BDTcut_wp97p5; - float bdtWP_XGB = 1. / (1. + std::sqrt((1. - bdtWP_MVA) / (1. + bdtWP_MVA))); // Convert WP value from ROOT.TMVA to XGboost + std::sort(candidates.begin(), + candidates.end(), + [](const CompositeCandidate &a, const CompositeCandidate &b) -> bool { return a.dpt < b.dpt; }); + unsigned int nCandPerCluster = std::min(candidates.size(), cfg.nCompCandPerCluster); + if (nCandPerCluster == 0) + continue; + float maxScore = -999; int ibest = -1; - for(unsigned int icand = 0; icand < nCandPerCluster; icand++) { + for (unsigned int icand = 0; icand < nCandPerCluster; icand++) { auto &cand = candidates[icand]; - std::vector emcalo_sel = emcalo; + const std::vector &emcalo_sel = emcalo; float score = compute_composite_score(cand, emcalo_sel, track, cfg.compIDparams); - if(score > maxScore) { - // if((score > bdtWP_XGB) && (score > maxScore)) { + if ((score > cfg.compIDparams.bdtScore_loose_wp) && (score > maxScore)) { maxScore = score; ibest = icand; } } - if(ibest != -1) { + if (ibest != -1) { emCalo2tk[ic] = candidates[ibest].track_idx; emCaloTkBdtScore[ic] = maxScore; } } } - float PFTkEGAlgoEmulator::compute_composite_score(CompositeCandidate &cand, const std::vector &emcalo, const std::vector &track, @@ -270,30 +233,28 @@ float PFTkEGAlgoEmulator::compute_composite_score(CompositeCandidate &cand, // Call and normalize input feature values, then cast to ap_fixed. // Note that for some features (e.g. track pT) we call the floating point representation, but that's already quantized! // Several other features, such as chi2 or most cluster features, are not quantized before casting them to ap_fixed. - ap_fixed<21,12,AP_RND_CONV,AP_SAT> hoe = calo.hwHoe; - ap_fixed<21,12,AP_RND_CONV,AP_SAT> tkpt = tk.hwPt; - ap_fixed<21,12,AP_RND_CONV,AP_SAT> srrtot = calo.hwSrrTot; - ap_fixed<21,12,AP_RND_CONV,AP_SAT> deta = tk.hwEta - calo.hwEta; - ap_fixed<18,9> calo_invPt = invert_with_shift, 1024>(calo.hwPt); // TODO: this is a guess - ap_fixed<21,12,AP_RND_CONV,AP_SAT> dpt = tk.hwPt * calo_invPt; - ap_fixed<21,12,AP_RND_CONV,AP_SAT> meanz = calo.hwMeanZ; - ap_fixed<21,12,AP_RND_CONV,AP_SAT> dphi = tk.hwPhi - calo.hwPhi; - ap_fixed<21,12,AP_RND_CONV,AP_SAT> chi2 = tk.hwChi2; - ap_fixed<21,12,AP_RND_CONV,AP_SAT> tkz0 = tk.hwZ0; - ap_fixed<21,12,AP_RND_CONV,AP_SAT> nstubs = tk.hwStubs; - + bdt_feature_t hoe = calo.hwHoe; + bdt_feature_t tkpt = tk.hwPt; + bdt_feature_t srrtot = calo.hwSrrTot; + bdt_feature_t deta = tk.hwEta - calo.hwEta; + ap_fixed<18, 9> calo_invPt = invert_with_shift, 1024>(calo.hwPt); // TODO: this is a guess + bdt_feature_t dpt = tk.hwPt * calo_invPt; + bdt_feature_t meanz = calo.hwMeanZ; + bdt_feature_t dphi = tk.hwPhi - calo.hwPhi; + bdt_feature_t chi2 = tk.hwRedChi2RPhi; + bdt_feature_t tkz0 = tk.hwZ0; + bdt_feature_t nstubs = tk.hwStubs; + // Run BDT inference - std::vector> inputs = { hoe, tkpt, srrtot, deta, dpt, meanz, dphi, chi2, tkz0, nstubs } ; - std::vector> bdt_score = composite_bdt_->decision_function(inputs); + std::vector inputs = {hoe, tkpt, srrtot, deta, dpt, meanz, dphi, chi2, tkz0, nstubs}; + std::vector> bdt_score = composite_bdt_->decision_function(inputs); float bdt_score_CON = bdt_score[0]; - float bdt_score_XGB = 1/(1+exp(-bdt_score_CON)); // Map Conifer score to XGboost score. (same as scipy.expit) + float bdt_score_XGB = 1 / (1 + exp(-bdt_score_CON)); // Map Conifer score to XGboost score. (same as scipy.expit) - // std::cout<<"BDT score of composite candidate = "< &emcalo, std::vector &emcalo_sel) const { @@ -319,7 +280,6 @@ void PFTkEGAlgoEmulator::run(const PFInputRegion &in, OutputRegion &out) const { << std::endl; } } - std::cout<<"running"< emCalo2tk(emcalo_sel.size(), -1); std::vector emCaloTkBdtScore(emcalo_sel.size(), -999); - std::cout<<"about to start matching"< egobjs; std::vector egeleobjs; @@ -475,9 +434,14 @@ EGIsoEleObjEmu &PFTkEGAlgoEmulator::addEGIsoEleToPF(std::vector egiso.hwPhi = calo.hwPhi; unsigned int egHwQual = hwQual; if (cfg.doEndcapHwQual) { - // 1. zero-suppress the loose EG-ID (bit 1) - // 2. for now use the standalone tight definition (bit 0) to set the tight point for eles (bit 1) - egHwQual = (hwQual & 0x9) | ((hwQual & 0x1) << 1); + if (cfg.doCompositeTkEle) { + // tight ele WP is set for tight BDT score + egHwQual = (hwQual & 0x9) | ((bdtScore >= cfg.compIDparams.bdtScore_tight_wp) << 1); + } else { + // 1. zero-suppress the loose EG-ID (bit 1) + // 2. for now use the standalone tight definition (bit 0) to set the tight point for eles (bit 1) + egHwQual = (hwQual & 0x9) | ((hwQual & 0x1) << 1); + } } egiso.hwQual = egHwQual; egiso.hwDEta = track.hwVtxEta() - egiso.hwEta; diff --git a/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp index 48d8fc85fb2cc..15bf8a01bca3f 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp @@ -13,7 +13,6 @@ l1ct::HadCaloObjEmu l1ct::HgcalClusterDecoderEmulator::decode(const ap_uint<256> // FIXME: we use a spare space in the word for hoe which is not in the current interface ap_uint<12> w_hoe = in(127, 116); - l1ct::HadCaloObjEmu out; out.clear(); out.hwPt = w_pt * l1ct::pt_t(l1ct::Scales::INTPT_LSB); @@ -24,8 +23,10 @@ l1ct::HadCaloObjEmu l1ct::HgcalClusterDecoderEmulator::decode(const ap_uint<256> // FIXME: variables use by composite-ID need to be added here out.hwSrrTot = w_srrtot * l1ct::srrtot_t(l1ct::Scales::SRRTOT_LSB); - out.hwMeanZ = l1ct::meanz_t(w_meanz - l1ct::Scales::MEANZ_SCALE); - calo.hwHoe = w_hoe * l1ct::hoe_t(l1ct::Scales::HOE_LSB); + out.hwMeanZ = (w_meanz == 0) ? l1ct::meanz_t(0) : l1ct::meanz_t(w_meanz - l1ct::meanz_t(l1ct::Scales::MEANZ_OFFSET)); + out.hwHoe = w_hoe * l1ct::hoe_t(l1ct::Scales::HOE_LSB); + + // std::cout << "[HadCaloObjEmu] meanz in: " << w_meanz << " out: " << out.hwMeanZ << std::endl; return out; } diff --git a/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/tkinput_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/tkinput_ref.cpp index 4f707e49cc3cc..ca61b38eb3eb6 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/tkinput_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/tkinput_ref.cpp @@ -37,7 +37,8 @@ namespace { l1ct::TrackInputEmulator::TrackInputEmulator(const edm::ParameterSet &iConfig) : TrackInputEmulator(parseRegion(iConfig.getParameter("region")), parseEncoding(iConfig.getParameter("trackWordEncoding")), - iConfig.getParameter("bitwiseAccurate")) { + iConfig.getParameter("bitwiseAccurate"), + iConfig.getParameter("slimDataFormat")) { if (region_ != Region::Endcap && region_ != Region::Barrel) { edm::LogError("TrackInputEmulator") << "region '" << iConfig.getParameter("region") << "' is not yet supported"; @@ -81,10 +82,11 @@ l1ct::TrackInputEmulator::TrackInputEmulator(const edm::ParameterSet &iConfig) #endif -l1ct::TrackInputEmulator::TrackInputEmulator(Region region, Encoding encoding, bool bitwise) +l1ct::TrackInputEmulator::TrackInputEmulator(Region region, Encoding encoding, bool bitwise, bool slim) : region_(region), encoding_(encoding), bitwise_(bitwise), + slim_(slim), rInvToPt_(31199.5), phiScale_(0.00038349520), z0Scale_(0.00999469), @@ -100,7 +102,8 @@ l1ct::TrackInputEmulator::TrackInputEmulator(Region region, Encoding encoding, b std::pair l1ct::TrackInputEmulator::decodeTrack(ap_uint<96> tkword, const l1ct::PFRegionEmu §or, - bool bitwise) const { + bool bitwise, + bool slim) const { l1ct::TkObjEmu ret; ret.clear(); auto z0 = signedZ0(tkword); @@ -162,7 +165,6 @@ std::pair l1ct::TrackInputEmulator::decodeTrack(ap_uint<96 fDEta = floatDEtaHGCal(z0, Rinv, tanl); fDPhi = floatDPhiHGCal(z0, Rinv, tanl); } - ret.hwDPhi = std::round(fDPhi); ret.hwDEta = std::round(fDEta); ret.hwPhi = std::round(fvtxPhi - fDPhi * ret.intCharge()); @@ -171,6 +173,14 @@ std::pair l1ct::TrackInputEmulator::decodeTrack(ap_uint<96 ret.hwZ0 = l1ct::Scales::makeZ0(floatZ0(z0)); } + if (!slim) { + ap_uint<7> w_hitPattern = tkword(15, 9); + ret.hwStubs = countSetBits(w_hitPattern); + ret.hwRedChi2RZ = tkword(35, 32); + ret.hwRedChi2RPhi = tkword(67, 64); + ret.hwRedChi2Bend = tkword(18, 16); + } + oksel = ret.hwQuality != 0; } return std::make_pair(ret, oksel); diff --git a/L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_dumpFiles_cfg.py b/L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_dumpFiles_cfg.py index edc30246e11c4..c8eddfb3b029c 100644 --- a/L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_dumpFiles_cfg.py +++ b/L1Trigger/Phase2L1ParticleFlow/test/make_l1ctLayer1_dumpFiles_cfg.py @@ -7,7 +7,7 @@ process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") process.load("FWCore.MessageLogger.MessageLogger_cfi") process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True), allowUnscheduled = cms.untracked.bool(False) ) -process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(100)) +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(1000)) process.MessageLogger.cerr.FwkReport.reportEvery = 1 process.source = cms.Source("PoolSource", From ecdb9d5f9d03ea5304ff01464cf5c43f282b856f Mon Sep 17 00:00:00 2001 From: Peter Eduard Meiring Date: Tue, 28 Mar 2023 10:53:40 +0200 Subject: [PATCH 09/13] new model with srrtot rescaled, binned chi2 variables and without z0 --- DataFormats/L1TParticleFlow/interface/datatypes.h | 2 +- .../Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/DataFormats/L1TParticleFlow/interface/datatypes.h b/DataFormats/L1TParticleFlow/interface/datatypes.h index 985e5ea6a4a9c..d94ad8ec4e0bc 100644 --- a/DataFormats/L1TParticleFlow/interface/datatypes.h +++ b/DataFormats/L1TParticleFlow/interface/datatypes.h @@ -207,7 +207,7 @@ namespace l1ct { inline iso_t makeIso(float iso) { return iso_t(0.25 * round(iso * 4)); } inline int makeDR2FromFloatDR(float dr) { return ceil(dr * dr / ETAPHI_LSB / ETAPHI_LSB); } - inline srrtot_t makeSrrTot(float var) { return srrtot_t(SRRTOT_LSB * round(var / SRRTOT_LSB)); }; + inline srrtot_t makeSrrTot(float var) { return srrtot_t(SRRTOT_LSB * round(var * pow(2,6) / SRRTOT_LSB)); }; inline meanz_t makeMeanZ(float var) { return round(var - MEANZ_OFFSET); }; inline hoe_t makeHoe(float var) { return hoe_t(HOE_LSB * round(var / HOE_LSB)); }; diff --git a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp index 67d43529efcdd..5109cb37f69fd 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp @@ -241,12 +241,13 @@ float PFTkEGAlgoEmulator::compute_composite_score(CompositeCandidate &cand, bdt_feature_t dpt = tk.hwPt * calo_invPt; bdt_feature_t meanz = calo.hwMeanZ; bdt_feature_t dphi = tk.hwPhi - calo.hwPhi; - bdt_feature_t chi2 = tk.hwRedChi2RPhi; - bdt_feature_t tkz0 = tk.hwZ0; bdt_feature_t nstubs = tk.hwStubs; + bdt_feature_t chi2rphi = tk.hwRedChi2RPhi; + bdt_feature_t chi2rz = tk.hwRedChi2RZ; + bdt_feature_t chi2bend = tk.hwRedChi2Bend; // Run BDT inference - std::vector inputs = {hoe, tkpt, srrtot, deta, dpt, meanz, dphi, chi2, tkz0, nstubs}; + std::vector inputs = {tkpt, hoe, srrtot, deta, dphi, dpt, meanz, nstubs, chi2rphi, chi2rz, chi2bend}; std::vector> bdt_score = composite_bdt_->decision_function(inputs); float bdt_score_CON = bdt_score[0]; From af1b5b24d5812a2f326d24cb7ac32f83784e6cdd Mon Sep 17 00:00:00 2001 From: Gianluca Date: Tue, 28 Mar 2023 14:31:32 +0200 Subject: [PATCH 10/13] Consolidate quantization of input features + cosmetics. --- .../L1TCorrelator/interface/TkElectron.h | 7 +++--- DataFormats/L1TCorrelator/src/classes_def.xml | 3 ++- .../L1TParticleFlow/interface/datatypes.h | 9 ++++---- .../interface/layer1_emulator.h | 4 ++-- .../L1TParticleFlow/interface/layer1_objs.h | 2 +- .../interface/common/inversion.h | 2 +- .../interface/egamma/pftkegalgo_ref.h | 1 + .../plugins/L1TCorrelatorLayer1Producer.cc | 7 +++--- .../plugins/L1TCtL2EgProducer.cc | 4 +++- .../python/l1TkEgAlgoEmulator_cfi.py | 11 ++++----- .../src/egamma/pftkegalgo_ref.cpp | 23 ++++++++----------- 11 files changed, 34 insertions(+), 39 deletions(-) diff --git a/DataFormats/L1TCorrelator/interface/TkElectron.h b/DataFormats/L1TCorrelator/interface/TkElectron.h index 9d66291cbc939..c5b70d79598fb 100644 --- a/DataFormats/L1TCorrelator/interface/TkElectron.h +++ b/DataFormats/L1TCorrelator/interface/TkElectron.h @@ -39,19 +39,18 @@ namespace l1t { float trkzVtx() const { return trkzVtx_; } double trackCurvature() const { return trackCurvature_; } - float compositeBdtScore() const { return compositeBdtScore_; } + float idScore() const { return idScore_; } // ---------- member functions --------------------------- void setTrkzVtx(float TrkzVtx) { trkzVtx_ = TrkzVtx; } void setTrackCurvature(double trackCurvature) { trackCurvature_ = trackCurvature; } - void setCompositeBdtScore(float score) { compositeBdtScore_ = score; } - + void setIdScore(float score) { idScore_ = score; } private: edm::Ptr trkPtr_; float trkzVtx_; double trackCurvature_; - float compositeBdtScore_; + float idScore_; }; } // namespace l1t #endif diff --git a/DataFormats/L1TCorrelator/src/classes_def.xml b/DataFormats/L1TCorrelator/src/classes_def.xml index 47267dc2b8eaa..723849343a5a9 100644 --- a/DataFormats/L1TCorrelator/src/classes_def.xml +++ b/DataFormats/L1TCorrelator/src/classes_def.xml @@ -44,7 +44,8 @@ - + + diff --git a/DataFormats/L1TParticleFlow/interface/datatypes.h b/DataFormats/L1TParticleFlow/interface/datatypes.h index d94ad8ec4e0bc..be802bc540b6e 100644 --- a/DataFormats/L1TParticleFlow/interface/datatypes.h +++ b/DataFormats/L1TParticleFlow/interface/datatypes.h @@ -156,8 +156,9 @@ namespace l1ct { constexpr float DXY_LSB = 0.05; constexpr float PUPPIW_LSB = 1.0 / 256; constexpr float MEANZ_OFFSET = 320.; - constexpr float SRRTOT_LSB = pow(2, -9); - constexpr float HOE_LSB = pow(2, -5); + constexpr float SRRTOT_LSB = 0.0019531250; // pow(2, -9) + constexpr unsigned int SRRTOT_SCALE = 64; // pow(2, 6) + constexpr float HOE_LSB = 0.031250000; // pow(2, -5) inline float floatPt(pt_t pt) { return pt.to_float(); } inline float floatPt(dpt_t pt) { return pt.to_float(); } @@ -174,7 +175,7 @@ namespace l1ct { inline float floatDxy(dxy_t dxy) { return dxy.to_float() * DXY_LSB; } inline float floatPuppiW(puppiWgt_t puppiw) { return puppiw.to_float() * PUPPIW_LSB; } inline float floatIso(iso_t iso) { return iso.to_float(); } - inline float floatSrrTot(srrtot_t srrtot) { return srrtot.to_float(); }; + inline float floatSrrTot(srrtot_t srrtot) { return srrtot.to_float() / SRRTOT_SCALE; }; inline float floatMeanZ(meanz_t meanz) { return meanz + MEANZ_OFFSET; }; inline float floatHoe(hoe_t hoe) { return hoe.to_float(); }; @@ -207,7 +208,7 @@ namespace l1ct { inline iso_t makeIso(float iso) { return iso_t(0.25 * round(iso * 4)); } inline int makeDR2FromFloatDR(float dr) { return ceil(dr * dr / ETAPHI_LSB / ETAPHI_LSB); } - inline srrtot_t makeSrrTot(float var) { return srrtot_t(SRRTOT_LSB * round(var * pow(2,6) / SRRTOT_LSB)); }; + inline srrtot_t makeSrrTot(float var) { return srrtot_t(SRRTOT_LSB * round(var * SRRTOT_SCALE / SRRTOT_LSB)); }; inline meanz_t makeMeanZ(float var) { return round(var - MEANZ_OFFSET); }; inline hoe_t makeHoe(float var) { return hoe_t(HOE_LSB * round(var / HOE_LSB)); }; diff --git a/DataFormats/L1TParticleFlow/interface/layer1_emulator.h b/DataFormats/L1TParticleFlow/interface/layer1_emulator.h index 944b2211e5f15..4a7709bbf52e0 100644 --- a/DataFormats/L1TParticleFlow/interface/layer1_emulator.h +++ b/DataFormats/L1TParticleFlow/interface/layer1_emulator.h @@ -195,7 +195,7 @@ namespace l1ct { const l1t::PFTrack *srcTrack = nullptr; // we use an index to the standalone object needed to retrieve a Ref when putting int sta_idx; - float bdtScore; + float idScore; bool read(std::fstream &from); bool write(std::fstream &to) const; void clear() { @@ -203,7 +203,7 @@ namespace l1ct { srcCluster = nullptr; srcTrack = nullptr; sta_idx = -1; - bdtScore = -999; + idScore = -999; clearIsoVars(); } diff --git a/DataFormats/L1TParticleFlow/interface/layer1_objs.h b/DataFormats/L1TParticleFlow/interface/layer1_objs.h index 28a008eba1207..48af69b993de1 100644 --- a/DataFormats/L1TParticleFlow/interface/layer1_objs.h +++ b/DataFormats/L1TParticleFlow/interface/layer1_objs.h @@ -175,7 +175,7 @@ namespace l1ct { stub_t hwStubs; redChi2Bin_t hwRedChi2RZ; // 4 bits redChi2Bin_t hwRedChi2RPhi; // 4 bits - //FIXME: is this actually filled? 3 bits would be enough + //FIXME: 3 bits would be enough redChi2Bin_t hwRedChi2Bend; // 4 bits enum TkQuality { PFLOOSE = 1, PFTIGHT = 2 }; diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h b/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h index ab91b6ff28fed..8628a4ca1af24 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h @@ -43,7 +43,7 @@ table_t invert_with_shift(in_t in) { // find the first '1' in the denominator int msb = 0; for (int b = 0; b < in.width; b++) { -// #pragma HLS unroll + // #pragma HLS unroll if (in[b]) msb = b; } diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h index 521722b9c5a19..363a4130f9f34 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h @@ -159,6 +159,7 @@ namespace l1ct { bool writeEgSta() const { return cfg.writeEgSta; } typedef ap_fixed<21, 12, AP_RND_CONV, AP_SAT> bdt_feature_t; + typedef ap_fixed<12, 3, AP_RND_CONV, AP_SAT> bdt_score_t; private: void link_emCalo2emCalo(const std::vector &emcalo, std::vector &emCalo2emCalo) const; diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc index 15d0a72e83b50..619ef495249a6 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc @@ -38,7 +38,6 @@ #include "DataFormats/L1Trigger/interface/EGamma.h" #include "DataFormats/L1TCorrelator/interface/TkEm.h" #include "DataFormats/L1TCorrelator/interface/TkEmFwd.h" -#include "DataFormats/L1THGCal/interface/HGCalMulticluster.h" //-------------------------------------------------------------------------------------------------- class L1TCorrelatorLayer1Producer : public edm::stream::EDProducer<> { @@ -698,8 +697,8 @@ void L1TCorrelatorLayer1Producer::addRawHgcalCluster(l1ct::DetectorSector w_eta = round(sec.region.localEta(c.eta()) / ETAPHI_LSB); ap_int<9> w_phi = round(sec.region.localPhi(c.phi()) / ETAPHI_LSB); ap_uint<10> w_qual = c.hwQual(); - - ap_uint<13> w_srrtot = round(c.sigmaRR() / l1ct::Scales::SRRTOT_LSB); + // FIXME: this is an arbitrary choice to keep the rounding consistent with the "addDecodedHadCalo" one + ap_uint<13> w_srrtot = round(c.sigmaRR() * l1ct::Scales::SRRTOT_SCALE / l1ct::Scales::SRRTOT_LSB); ap_uint<12> w_meanz = round(c.absZBarycenter()); // FIXME: the calibration can actually make hoe become negative....we add a small protection for now // We use ap_ufixed to handle saturation and rounding @@ -1101,7 +1100,7 @@ void L1TCorrelatorLayer1Producer::putEgObjects(edm::Event &iEvent, tkele.setHwQual(egele.hwQual); tkele.setPFIsol(egele.floatRelIso(l1ct::EGIsoEleObjEmu::IsoType::PfIso)); tkele.setEgBinaryWord(egele.pack()); - tkele.setCompositeBdtScore(egele.bdtScore); + tkele.setIdScore(egele.idScore); tkeles->push_back(tkele); nele_obj.push_back(tkeles->size() - 1); } diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCtL2EgProducer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCtL2EgProducer.cc index 99b0251eff85b..54399d8a497a3 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCtL2EgProducer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCtL2EgProducer.cc @@ -351,7 +351,8 @@ void L1TCtL2EgProducer::convertToEmu(const l1t::TkElectron &tkele, emu.setHwIso(EGIsoEleObjEmu::IsoType::PfIso, l1ct::Scales::makeIso(tkele.pfIsol() * tkele.pt())); emu.setHwIso(EGIsoEleObjEmu::IsoType::PuppiIso, l1ct::Scales::makeIso(tkele.puppiIsol() * tkele.pt())); // std::cout << "[convertToEmu] TkEle pt: " << emu.hwPt << " eta: " << emu.hwEta << " phi: " << emu.hwPhi << " staidx: " << emu.sta_idx << std::endl; - + // FIXME: this is temporary while waiting to move the BDT score to the FW object + emu.idScore = tkele.idScore(); boarOut.egelectron.push_back(emu); } @@ -421,6 +422,7 @@ l1t::TkElectron L1TCtL2EgProducer::convertFromEmu(const l1ct::EGIsoEleObjEmu &eg tkele.setPFIsol(egele.floatRelIso(l1ct::EGIsoEleObjEmu::IsoType::PfIso)); tkele.setPuppiIsol(egele.floatRelIso(l1ct::EGIsoEleObjEmu::IsoType::PuppiIso)); tkele.setEgBinaryWord(gteg.pack()); + tkele.setIdScore(egele.idScore); return tkele; } diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py index 96693cc0f86ae..7f7661a03d588 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py @@ -1,8 +1,4 @@ import FWCore.ParameterSet.Config as cms -import math - -def BDTScore_tmva2xgboost(score): - return 1. / (1. + math.sqrt((1. - score) / (1. + score))) tkEgAlgoParameters = cms.PSet( nTRACK=cms.uint32(50), # very large numbers for first test @@ -57,9 +53,10 @@ def BDTScore_tmva2xgboost(score): doCompositeTkEle=cms.bool(False), nCompCandPerCluster=cms.uint32(3), compositeParametersTkEle=cms.PSet( - bdt_loose_wp=cms.double(0.), - bdt_tight_wp=cms.double(BDTScore_tmva2xgboost(0.9826955)), - conifer_model=cms.string("L1Trigger/Phase2L1ParticleFlow/data/compositeID.json") + # the working points are cuts on BDT output logits log(p/1-p) + loose_wp=cms.double(-4.), + tight_wp=cms.double(-1.029296875), + model=cms.string("L1Trigger/Phase2L1ParticleFlow/data/compositeID.json") ), ) diff --git a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp index 5109cb37f69fd..3f78f3529253f 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp @@ -52,9 +52,9 @@ l1ct::PFTkEGAlgoEmuConfig::IsoParameters::IsoParameters(const edm::ParameterSet pset.getParameter("dRMax")) {} l1ct::PFTkEGAlgoEmuConfig::CompIDParameters::CompIDParameters(const edm::ParameterSet &pset) - : CompIDParameters(pset.getParameter("bdt_loose_wp"), - pset.getParameter("bdt_tight_wp"), - pset.getParameter("conifer_model")) {} + : CompIDParameters(pset.getParameter("loose_wp"), + pset.getParameter("tight_wp"), + pset.getParameter("model")) {} #endif @@ -66,7 +66,7 @@ PFTkEGAlgoEmulator::PFTkEGAlgoEmulator(const PFTkEGAlgoEmuConfig &config) #else auto resolvedFileName = cfg.compIDparams.conifer_model; #endif - composite_bdt_ = new conifer::BDT, false>(resolvedFileName); + composite_bdt_ = new conifer::BDT(resolvedFileName); } } @@ -230,14 +230,12 @@ float PFTkEGAlgoEmulator::compute_composite_score(CompositeCandidate &cand, const auto &calo = emcalo[cand.cluster_idx]; const auto &tk = track[cand.track_idx]; - // Call and normalize input feature values, then cast to ap_fixed. - // Note that for some features (e.g. track pT) we call the floating point representation, but that's already quantized! - // Several other features, such as chi2 or most cluster features, are not quantized before casting them to ap_fixed. + // Prepare the input features bdt_feature_t hoe = calo.hwHoe; bdt_feature_t tkpt = tk.hwPt; bdt_feature_t srrtot = calo.hwSrrTot; bdt_feature_t deta = tk.hwEta - calo.hwEta; - ap_fixed<18, 9> calo_invPt = invert_with_shift, 1024>(calo.hwPt); // TODO: this is a guess + ap_ufixed<18, 0> calo_invPt = invert_with_shift, 1024>(calo.hwPt); bdt_feature_t dpt = tk.hwPt * calo_invPt; bdt_feature_t meanz = calo.hwMeanZ; bdt_feature_t dphi = tk.hwPhi - calo.hwPhi; @@ -248,12 +246,9 @@ float PFTkEGAlgoEmulator::compute_composite_score(CompositeCandidate &cand, // Run BDT inference std::vector inputs = {tkpt, hoe, srrtot, deta, dphi, dpt, meanz, nstubs, chi2rphi, chi2rz, chi2bend}; - std::vector> bdt_score = composite_bdt_->decision_function(inputs); + std::vector bdt_score = composite_bdt_->decision_function(inputs); - float bdt_score_CON = bdt_score[0]; - float bdt_score_XGB = 1 / (1 + exp(-bdt_score_CON)); // Map Conifer score to XGboost score. (same as scipy.expit) - - return bdt_score_XGB; + return bdt_score[0]; } void PFTkEGAlgoEmulator::sel_emCalo(unsigned int nmax_sel, @@ -451,7 +446,7 @@ EGIsoEleObjEmu &PFTkEGAlgoEmulator::addEGIsoEleToPF(std::vector egiso.hwCharge = track.hwCharge; egiso.srcCluster = calo.src; egiso.srcTrack = track.src; - egiso.bdtScore = bdtScore; + egiso.idScore = bdtScore; egobjs.push_back(egiso); if (debug_ > 2) From 11b293f1701f5d87dfbdda352d6353697c6661b0 Mon Sep 17 00:00:00 2001 From: Gianluca Date: Mon, 3 Apr 2023 11:34:18 +0200 Subject: [PATCH 11/13] Update tight WP to use 90%eff score. --- .../interface/common/inversion.h | 109 +++++++++--------- .../interface/egamma/pftkegalgo_ref.h | 6 +- .../plugins/L1TCorrelatorLayer1Producer.cc | 4 +- .../python/l1TkEgAlgoEmulator_cfi.py | 2 +- .../python/l1ctLayer1_cff.py | 55 ++++++++- .../python/l1ctLayer2EG_cff.py | 40 ++++++- .../src/egamma/pftkegalgo_ref.cpp | 2 +- .../src/l1-converters/hgcalinputt_ref.cpp | 3 - 8 files changed, 154 insertions(+), 67 deletions(-) diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h b/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h index 8628a4ca1af24..bf784f7bdbb5a 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h @@ -1,60 +1,63 @@ #ifndef CC_INVERSION_H__ #define CC_INVERSION_H__ -constexpr int ceillog2(int x) { return (x <= 2) ? 1 : 1 + ceillog2((x + 1) / 2); } - -template -inline float real_val_from_idx(unsigned i) { - // Treat the index as the top N bits - static constexpr int NB = ceillog2(N); // number of address bits for table - data_T x(0); - // The MSB of 1 is implicit in the table - x[x.width - 1] = 1; - // So we can use the next NB bits for real data - x(x.width - 2, x.width - NB - 1) = i; - return (float)x; -} - -template -inline unsigned idx_from_real_val(data_T x) { - // Slice the top N bits to get an index into the table - static constexpr int NB = ceillog2(N); // number of address bits for table - // Slice the top-1 NB bits of the value - // the MSB of '1' is implicit, so only slice below that - ap_uint y = x(x.width - 2, x.width - NB - 1); - return (unsigned)y(NB - 1, 0); -} - -template -void init_invert_table(table_T table_out[N]) { - // The template data_T is the data type used to address the table - for (unsigned i = 0; i < N; i++) { - float x = real_val_from_idx(i); - table_T inv_x = 1 / x; - table_out[i] = inv_x; +namespace l1ct { + + constexpr int ceillog2(int x) { return (x <= 2) ? 1 : 1 + ceillog2((x + 1) / 2); } + + template + inline float real_val_from_idx(unsigned i) { + // Treat the index as the top N bits + static constexpr int NB = ceillog2(N); // number of address bits for table + data_T x(0); + // The MSB of 1 is implicit in the table + x[x.width - 1] = 1; + // So we can use the next NB bits for real data + x(x.width - 2, x.width - NB - 1) = i; + return (float)x; + } + + template + inline unsigned idx_from_real_val(data_T x) { + // Slice the top N bits to get an index into the table + static constexpr int NB = ceillog2(N); // number of address bits for table + // Slice the top-1 NB bits of the value + // the MSB of '1' is implicit, so only slice below that + ap_uint y = x(x.width - 2, x.width - NB - 1); + return (unsigned)y(NB - 1, 0); } -} - -template -table_t invert_with_shift(in_t in) { - table_t inv_table[N]; - init_invert_table(inv_table); - - // find the first '1' in the denominator - int msb = 0; - for (int b = 0; b < in.width; b++) { - // #pragma HLS unroll - if (in[b]) - msb = b; + + template + void init_invert_table(table_T table_out[N]) { + // The template data_T is the data type used to address the table + for (unsigned i = 0; i < N; i++) { + float x = real_val_from_idx(i); + table_T inv_x = 1 / x; + table_out[i] = inv_x; + } + } + + template + table_t invert_with_shift(in_t in) { + table_t inv_table[N]; + init_invert_table(inv_table); + + // find the first '1' in the denominator + int msb = 0; + for (int b = 0; b < in.width; b++) { + // #pragma HLS unroll + if (in[b]) + msb = b; + } + // shift up the denominator such that the left-most bit (msb) is '1' + in_t in_shifted = in << (in.width - msb - 1); + // lookup the inverse of the shifted input + int idx = idx_from_real_val(in_shifted); + table_t inv_in = inv_table[idx]; + // shift the output back + table_t out = inv_in << (in.width - msb - 1); + return out; } - // shift up the denominator such that the left-most bit (msb) is '1' - in_t in_shifted = in << (in.width - msb - 1); - // lookup the inverse of the shifted input - int idx = idx_from_real_val(in_shifted); - table_t inv_in = inv_table[idx]; - // shift the output back - table_t out = inv_in << (in.width - msb - 1); - return out; -} +} // namespace l1ct #endif \ No newline at end of file diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h index 363a4130f9f34..560325e46cb0c 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h @@ -99,8 +99,7 @@ namespace l1ct { bool doPfIso = false, EGIsoEleObjEmu::IsoType hwIsoTypeTkEle = EGIsoEleObjEmu::IsoType::TkIso, EGIsoObjEmu::IsoType hwIsoTypeTkEm = EGIsoObjEmu::IsoType::TkIsoPV, - // FIXME: maybe we round these? - const CompIDParameters &compIDparams = {0.7927004, 0.9826955, "compositeID.json"}, + const CompIDParameters &compIDparams = {-4, 0.214844, "compositeID.json"}, int debug = 0) : nTRACK(nTrack), @@ -115,9 +114,6 @@ namespace l1ct { emClusterPtMin(emClusterPtMin), dEtaMaxBrem(dEtaMaxBrem), dPhiMaxBrem(dPhiMaxBrem), - //absEtaBoundaries(std::move(absEtaBoundaries)), - //dEtaValues(std::move(dEtaValues)), - //dPhiValues(std::move(dPhiValues)), absEtaBoundaries(absEtaBoundaries), dEtaValues(dEtaValues), dPhiValues(dPhiValues), diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc index 619ef495249a6..f53e14e2a5fc7 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc @@ -697,10 +697,10 @@ void L1TCorrelatorLayer1Producer::addRawHgcalCluster(l1ct::DetectorSector w_eta = round(sec.region.localEta(c.eta()) / ETAPHI_LSB); ap_int<9> w_phi = round(sec.region.localPhi(c.phi()) / ETAPHI_LSB); ap_uint<10> w_qual = c.hwQual(); - // FIXME: this is an arbitrary choice to keep the rounding consistent with the "addDecodedHadCalo" one + // NOTE: this is an arbitrary choice to keep the rounding consistent with the "addDecodedHadCalo" one ap_uint<13> w_srrtot = round(c.sigmaRR() * l1ct::Scales::SRRTOT_SCALE / l1ct::Scales::SRRTOT_LSB); ap_uint<12> w_meanz = round(c.absZBarycenter()); - // FIXME: the calibration can actually make hoe become negative....we add a small protection for now + // NOTE: the calibration can actually make hoe become negative....we add a small protection for now // We use ap_ufixed to handle saturation and rounding ap_ufixed<10, 5, AP_RND_CONV, AP_SAT> w_hoe = c.hOverE(); diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py index 7f7661a03d588..a060853085cbd 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py @@ -55,7 +55,7 @@ compositeParametersTkEle=cms.PSet( # the working points are cuts on BDT output logits log(p/1-p) loose_wp=cms.double(-4.), - tight_wp=cms.double(-1.029296875), + tight_wp=cms.double(0.214844), model=cms.string("L1Trigger/Phase2L1ParticleFlow/data/compositeID.json") ), ) diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py index c1e1acd8cff18..ac92502e74128 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py @@ -294,8 +294,15 @@ writeRawHgcalCluster = cms.untracked.bool(True) ) + l1tLayer1HGCalExtended = l1tLayer1HGCal.clone(tracks = ('l1tPFTracksFromL1TracksExtended')) +l1tLayer1HGCalElliptic = l1tLayer1HGCal.clone( + tkEgAlgoParameters = l1tLayer1HGCal.tkEgAlgoParameters.clone( + doCompositeTkEle = False, + trkQualityPtMin = 10.) +) + l1tLayer1HGCalNoTK = cms.EDProducer("L1TCorrelatorLayer1Producer", tracks = cms.InputTag(''), muons = cms.InputTag('l1tSAMuonsGmt','promptSAMuons'), @@ -526,6 +533,50 @@ ) ) +l1tLayer1EGElliptic = cms.EDProducer( + "L1TEGMultiMerger", + tkElectrons=cms.VPSet( + cms.PSet( + instance=cms.string("L1TkEleEE"), + pfProducers=cms.VInputTag( + cms.InputTag("l1tLayer1HGCalElliptic", 'L1TkEle') + ) + ), + cms.PSet( + instance=cms.string("L1TkEleEB"), + pfProducers=cms.VInputTag( + cms.InputTag("l1tLayer1Barrel", 'L1TkEle') + ) + ) + ), + tkEms=cms.VPSet( + cms.PSet( + instance=cms.string("L1TkEmEE"), + pfProducers=cms.VInputTag( + cms.InputTag("l1tLayer1HGCalElliptic", 'L1TkEm'), + cms.InputTag("l1tLayer1HGCalNoTK", 'L1TkEm') + ) + ), + cms.PSet( + instance=cms.string("L1TkEmEB"), + pfProducers=cms.VInputTag( + cms.InputTag("l1tLayer1Barrel", 'L1TkEm') + ) + ) + ), + tkEgs=cms.VPSet( + cms.PSet( + instance=cms.string("L1EgEE"), + pfProducers=cms.VInputTag( + cms.InputTag("l1tLayer1HGCalElliptic", 'L1Eg'), + cms.InputTag("l1tLayer1HGCalNoTK", 'L1Eg') + ) + ) + ) +) + + + L1TLayer1TaskInputsTask = cms.Task( l1tPFClustersFromL1EGClusters, l1tPFClustersFromCombinedCaloHCal, @@ -544,5 +595,7 @@ l1tLayer1HF, l1tLayer1, l1tLayer1Extended, - l1tLayer1EG + l1tLayer1HGCalElliptic, + l1tLayer1EG, + l1tLayer1EGElliptic ) diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer2EG_cff.py b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer2EG_cff.py index f9d9196a1611e..c0e7943579dfa 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer2EG_cff.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer2EG_cff.py @@ -150,8 +150,46 @@ # ) ) +l1tLayer2EGElliptic = l1tLayer2EG.clone( + tkElectrons=cms.VPSet( + cms.PSet( + pfProducer=cms.InputTag("l1tLayer1HGCalElliptic", 'L1TkElePerBoard'), + channels=cms.vint32(3, 4) + ), + cms.PSet( + pfProducer=cms.InputTag("l1tLayer1Barrel", 'L1TkElePerBoard'), + channels=cms.vint32(0, 1, 2) + ), + ), + tkEms=cms.VPSet( + cms.PSet( + pfProducer=cms.InputTag("l1tLayer1HGCalElliptic", 'L1TkEmPerBoard'), + channels=cms.vint32(3, 4) + ), + cms.PSet( + pfProducer=cms.InputTag("l1tLayer1HGCalNoTK", 'L1TkEmPerBoard'), + channels=cms.vint32(-1) + ), + cms.PSet( + pfProducer=cms.InputTag("l1tLayer1Barrel", 'L1TkEmPerBoard'), + channels=cms.vint32(0, 1, 2) + ), + ), + tkEgs=cms.VPSet( + cms.PSet( + pfProducer=cms.InputTag("l1tLayer1HGCalElliptic", 'L1Eg'), + channels=cms.vint32(-1) + ), + cms.PSet( + pfProducer=cms.InputTag("l1tLayer1HGCalNoTK", 'L1Eg'), + channels=cms.vint32(-1) + ), + ), +) + L1TLayer2EGTask = cms.Task( l1tLayer2Deregionizer, - l1tLayer2EG + l1tLayer2EG, + l1tLayer2EGElliptic ) diff --git a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp index 3f78f3529253f..c149f57006096 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp @@ -235,7 +235,7 @@ float PFTkEGAlgoEmulator::compute_composite_score(CompositeCandidate &cand, bdt_feature_t tkpt = tk.hwPt; bdt_feature_t srrtot = calo.hwSrrTot; bdt_feature_t deta = tk.hwEta - calo.hwEta; - ap_ufixed<18, 0> calo_invPt = invert_with_shift, 1024>(calo.hwPt); + ap_ufixed<18, 0> calo_invPt = l1ct::invert_with_shift, 1024>(calo.hwPt); bdt_feature_t dpt = tk.hwPt * calo_invPt; bdt_feature_t meanz = calo.hwMeanZ; bdt_feature_t dphi = tk.hwPhi - calo.hwPhi; diff --git a/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp index 15bf8a01bca3f..4851eccc06b71 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp @@ -21,12 +21,9 @@ l1ct::HadCaloObjEmu l1ct::HgcalClusterDecoderEmulator::decode(const ap_uint<256> out.hwEmPt = w_empt * l1ct::pt_t(l1ct::Scales::INTPT_LSB); out.hwEmID = w_qual; - // FIXME: variables use by composite-ID need to be added here out.hwSrrTot = w_srrtot * l1ct::srrtot_t(l1ct::Scales::SRRTOT_LSB); out.hwMeanZ = (w_meanz == 0) ? l1ct::meanz_t(0) : l1ct::meanz_t(w_meanz - l1ct::meanz_t(l1ct::Scales::MEANZ_OFFSET)); out.hwHoe = w_hoe * l1ct::hoe_t(l1ct::Scales::HOE_LSB); - // std::cout << "[HadCaloObjEmu] meanz in: " << w_meanz << " out: " << out.hwMeanZ << std::endl; - return out; } From 454105dd138883e66addb9d1d79b5ce938af4015 Mon Sep 17 00:00:00 2001 From: Gianluca Date: Thu, 13 Apr 2023 16:41:23 +0200 Subject: [PATCH 12/13] use the 0.955 score for the loose WP --- L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py index a060853085cbd..4b267ff987db1 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py @@ -54,7 +54,7 @@ nCompCandPerCluster=cms.uint32(3), compositeParametersTkEle=cms.PSet( # the working points are cuts on BDT output logits log(p/1-p) - loose_wp=cms.double(-4.), + loose_wp=cms.double(-0.732422), tight_wp=cms.double(0.214844), model=cms.string("L1Trigger/Phase2L1ParticleFlow/data/compositeID.json") ), From efa052cd7afa7f5789bfb061936e7aaedf2054b7 Mon Sep 17 00:00:00 2001 From: Gianluca Date: Tue, 16 May 2023 18:18:24 +0200 Subject: [PATCH 13/13] implement review comments --- DataFormats/L1TCorrelator/src/classes_def.xml | 5 ++--- L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp | 4 ++-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/DataFormats/L1TCorrelator/src/classes_def.xml b/DataFormats/L1TCorrelator/src/classes_def.xml index 723849343a5a9..6b5acaf05cbc1 100644 --- a/DataFormats/L1TCorrelator/src/classes_def.xml +++ b/DataFormats/L1TCorrelator/src/classes_def.xml @@ -44,9 +44,8 @@ - - - + + diff --git a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp index c149f57006096..984dca69f36a5 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp @@ -185,9 +185,9 @@ void PFTkEGAlgoEmulator::link_emCalo2tk_composite(const PFRegionEmu &r, float d_phi = deltaPhi(tk.floatPhi(), calo.floatPhi()); float d_eta = tk.floatEta() - calo.floatEta(); // We only use it squared - float dR = std::sqrt((d_phi * d_phi) + (d_eta * d_eta)); + float dR2 = (d_phi * d_phi) + (d_eta * d_eta); - if (dR < 0.2) { + if (dR2 < 0.04) { // Only store indices, dR and dpT for now. The other quantities are computed only for the best nCandPerCluster. CompositeCandidate cand; cand.cluster_idx = ic;