diff --git a/DataFormats/L1TCorrelator/interface/TkElectron.h b/DataFormats/L1TCorrelator/interface/TkElectron.h index 9f1ac1eaeeff0..c5b70d79598fb 100644 --- a/DataFormats/L1TCorrelator/interface/TkElectron.h +++ b/DataFormats/L1TCorrelator/interface/TkElectron.h @@ -39,16 +39,18 @@ namespace l1t { float trkzVtx() const { return trkzVtx_; } double trackCurvature() const { return trackCurvature_; } - + float idScore() const { return idScore_; } // ---------- member functions --------------------------- void setTrkzVtx(float TrkzVtx) { trkzVtx_ = TrkzVtx; } void setTrackCurvature(double trackCurvature) { trackCurvature_ = trackCurvature; } + void setIdScore(float score) { idScore_ = score; } private: edm::Ptr trkPtr_; float trkzVtx_; double trackCurvature_; + float idScore_; }; } // namespace l1t #endif diff --git a/DataFormats/L1TCorrelator/src/classes_def.xml b/DataFormats/L1TCorrelator/src/classes_def.xml index 0685c6bc1356e..6b5acaf05cbc1 100644 --- a/DataFormats/L1TCorrelator/src/classes_def.xml +++ b/DataFormats/L1TCorrelator/src/classes_def.xml @@ -44,7 +44,8 @@ - + + diff --git a/DataFormats/L1TMuonPhase2/interface/Constants.h b/DataFormats/L1TMuonPhase2/interface/Constants.h index d3f80a3421948..ea8c2f4c0c5cb 100644 --- a/DataFormats/L1TMuonPhase2/interface/Constants.h +++ b/DataFormats/L1TMuonPhase2/interface/Constants.h @@ -54,7 +54,7 @@ namespace Phase2L1GMT { // Bitwidth for standalone muons to CL1 and GT const int BITSSAZ0 = 5; const int BITSSAD0 = 7; - const int BITSSAQUALITY = 4; + const int BITSSAQUAL = 4; // Bitwidth for dataformat to GT const int BITSGTPT = 16; @@ -62,14 +62,24 @@ namespace Phase2L1GMT { const int BITSGTETA = 14; const int BITSGTZ0 = 10; const int BITSGTD0 = 10; - const int BITSGTQUALITY = 8; + const int BITSGTQUAL = 8; const int BITSGTISO = 4; + const int BITSGTBETA = 4; + + // Bitwidth for Tau->3mu object + const int BITSTMPT = 8; + const int BITSTMPHI = 8; + const int BITSTMETA = 8; + const int BITSTMMASS2 = 8; + const int BITSTMTYPE = 6; + const int BITSTMIDX = 4; + const int BITSTMQUAL = 4; const float maxCurv_ = 0.00855; // 2 GeV pT Rinv is in cm const float maxPhi_ = 1.026; // relative to the center of the sector const float maxTanl_ = 8.0; - const float maxZ0_ = 30.; - const float maxD0_ = 15.4; + const float maxZ0_ = 25.6; + const float maxD0_ = 15.36; // Updated barrelLimit according to Karol, https://indico.cern.ch/event/1113802/#1-phase2-gmt-performance-and-i const int barrelLimit0_ = 1.4 / 0.00076699039 / 8; const int barrelLimit1_ = 1.1 / 0.00076699039 / 8; @@ -81,12 +91,32 @@ namespace Phase2L1GMT { const float LSBpt = 0.03125; const float LSBphi = 2. * M_PI / pow(2, BITSPHI); const float LSBeta = 2. * M_PI / pow(2, BITSETA); - const float LSBGTz0 = 2. * maxZ0_ / pow(2, BITSZ0); - const float LSBGTd0 = 2. * maxD0_ / pow(2, BITSD0); - const float LSBSAz0 = 1.875; - const float LSBSAd0 = 3.85; + const float LSBGTz0 = 0.05; // 0.5mm, in sync with GTT and Correlator + const float LSBGTd0 = 0.03; // from GT interface doc + const float LSBSAz0 = 1.6; // 0.05 * 32 cm, with range +- 25.6 + const float LSBSAd0 = 3.84; // 0.03 * 128 cm, with range +- 245.76 typedef ap_uint<64> wordtype; + typedef ap_uint<1> valid_gt_t; //valid + typedef ap_uint<1> q_gt_t; //charge + typedef ap_uint pt_gt_t; //pt of tracker muon + typedef ap_int phi_gt_t; //phi of tracker muon + typedef ap_int eta_gt_t; //eta of tracker muon + typedef ap_int z0_gt_t; //z0 of tracker muon + typedef ap_int d0_gt_t; //d0 of tracker muon + typedef ap_uint iso_gt_t; //isolation of tracker muon + typedef ap_uint beta_gt_t; //beta of tracker muon + typedef ap_uint qual_gt_t; //quality of tracker muon + + //Standalone muon datatype + typedef ap_uint<1> valid_sa_t; //valid + typedef ap_uint pt_sa_t; //pt of standalone muon + typedef ap_int phi_sa_t; //phi of standalone muon + typedef ap_int eta_sa_t; //eta of standalone muon + typedef ap_int z0_sa_t; //z0 of standalone muon + typedef ap_int d0_sa_t; //d0 of standalone muon + typedef ap_uint<1> q_sa_t; //charge of standalone muon + typedef ap_uint qual_sa_t; //quality of standalone muon inline uint64_t twos_complement(long long int v, uint bits) { uint64_t mask = (1 << bits) - 1; diff --git a/DataFormats/L1TMuonPhase2/interface/SAMuon.h b/DataFormats/L1TMuonPhase2/interface/SAMuon.h index d46d6cab3610c..12ea802dae04f 100644 --- a/DataFormats/L1TMuonPhase2/interface/SAMuon.h +++ b/DataFormats/L1TMuonPhase2/interface/SAMuon.h @@ -31,6 +31,16 @@ namespace l1t { const uint hwBeta() const { return hwBeta_; } void setBeta(uint beta) { hwBeta_ = beta; } + // For GT, returning ap_ type + const Phase2L1GMT::valid_sa_t apValid() const { return Phase2L1GMT::valid_sa_t(hwPt() > 0); }; + const Phase2L1GMT::pt_sa_t apPt() const { return Phase2L1GMT::pt_sa_t(hwPt()); }; + const Phase2L1GMT::phi_sa_t apPhi() const { return Phase2L1GMT::phi_sa_t(hwPhi()); }; + const Phase2L1GMT::eta_sa_t apEta() const { return Phase2L1GMT::eta_sa_t(hwEta()); }; + const Phase2L1GMT::z0_sa_t apZ0() const { return Phase2L1GMT::z0_sa_t(hwZ0()); }; + const Phase2L1GMT::d0_sa_t apD0() const { return Phase2L1GMT::d0_sa_t(hwD0()); }; + const Phase2L1GMT::q_sa_t apCharge() const { return Phase2L1GMT::q_sa_t(hwCharge()); }; + const Phase2L1GMT::qual_sa_t apQual() const { return Phase2L1GMT::qual_sa_t(hwQual()); }; + // For HLT const double phZ0() const { return Phase2L1GMT::LSBSAz0 * hwZ0(); } const double phD0() const { return Phase2L1GMT::LSBSAd0 * hwD0(); } diff --git a/DataFormats/L1TMuonPhase2/interface/TrackerMuon.h b/DataFormats/L1TMuonPhase2/interface/TrackerMuon.h index 78bb288685d37..1dd8cc90bf344 100644 --- a/DataFormats/L1TMuonPhase2/interface/TrackerMuon.h +++ b/DataFormats/L1TMuonPhase2/interface/TrackerMuon.h @@ -32,7 +32,7 @@ namespace l1t { ~TrackerMuon() override; const edm::Ptr& trkPtr() const { return trkPtr_; } - const edm::Ref& muonRef() const { return muRef_; } + const std::vector& muonRef() const { return muRef_; } const bool hwCharge() const { return hwCharge_; } const int hwZ0() const { return hwZ0_; } @@ -41,10 +41,22 @@ namespace l1t { const int hwIsoSumAp() const { return hwIsoSumAp_; } const uint hwBeta() const { return hwBeta_; } void setBeta(uint beta) { hwBeta_ = beta; } - void setMuonRef(const edm::Ref& p) { muRef_ = p; } + void setMuonRef(const std::vector& p) { muRef_ = p; } void setHwIsoSum(int isoSum) { hwIsoSum_ = isoSum; } void setHwIsoSumAp(int isoSum) { hwIsoSumAp_ = isoSum; } + // For GT, returning ap_ type + const Phase2L1GMT::valid_gt_t apValid() const { return Phase2L1GMT::valid_gt_t(hwPt() > 0); }; + const Phase2L1GMT::pt_gt_t apPt() const { return Phase2L1GMT::pt_gt_t(hwPt()); }; + const Phase2L1GMT::phi_gt_t apPhi() const { return Phase2L1GMT::phi_gt_t(hwPhi()); }; + const Phase2L1GMT::eta_gt_t apEta() const { return Phase2L1GMT::eta_gt_t(hwEta()); }; + const Phase2L1GMT::z0_gt_t apZ0() const { return Phase2L1GMT::z0_gt_t(hwZ0()); }; + const Phase2L1GMT::d0_gt_t apD0() const { return Phase2L1GMT::d0_gt_t(hwD0()); }; + const Phase2L1GMT::q_gt_t apCharge() const { return Phase2L1GMT::q_gt_t(hwCharge()); }; + const Phase2L1GMT::qual_gt_t apQual() const { return Phase2L1GMT::qual_gt_t(hwQual()); }; + const Phase2L1GMT::iso_gt_t apIso() const { return Phase2L1GMT::iso_gt_t(hwIso()); }; + const Phase2L1GMT::beta_gt_t apBeta() const { return Phase2L1GMT::beta_gt_t(hwBeta()); }; + // For HLT const double phZ0() const { return Phase2L1GMT::LSBGTz0 * hwZ0(); } const double phD0() const { return Phase2L1GMT::LSBGTd0 * hwD0(); } @@ -76,7 +88,7 @@ namespace l1t { //Store the eneryg sum for isolation with ap_type int hwIsoSumAp_; - edm::Ref muRef_; + std::vector muRef_; MuonStubRefVector stubs_; }; } // namespace l1t diff --git a/DataFormats/L1TMuonPhase2/src/classes_def.xml b/DataFormats/L1TMuonPhase2/src/classes_def.xml index b20feacd4c091..b3cbf00fcf430 100644 --- a/DataFormats/L1TMuonPhase2/src/classes_def.xml +++ b/DataFormats/L1TMuonPhase2/src/classes_def.xml @@ -8,7 +8,8 @@ - + + diff --git a/DataFormats/L1TParticleFlow/interface/PFCandidate.h b/DataFormats/L1TParticleFlow/interface/PFCandidate.h index c6d52246a6fa5..567eb6f80b7fb 100644 --- a/DataFormats/L1TParticleFlow/interface/PFCandidate.h +++ b/DataFormats/L1TParticleFlow/interface/PFCandidate.h @@ -48,9 +48,13 @@ 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 +74,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/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/PFJet.h b/DataFormats/L1TParticleFlow/interface/PFJet.h index e1cfc282659c0..5dd73fdc0bdec 100644 --- a/DataFormats/L1TParticleFlow/interface/PFJet.h +++ b/DataFormats/L1TParticleFlow/interface/PFJet.h @@ -5,6 +5,8 @@ #include "DataFormats/L1Trigger/interface/L1Candidate.h" #include "DataFormats/L1TParticleFlow/interface/PFCandidate.h" #include "DataFormats/Common/interface/Ptr.h" +#include "DataFormats/L1TParticleFlow/interface/jets.h" +#include "DataFormats/L1TParticleFlow/interface/gt_datatypes.h" namespace l1t { @@ -38,13 +40,24 @@ namespace l1t { edm::Ptr daughterPtr(size_type i) const { return constituents_[i]; } // Get and set the encodedJet_ bits. The Jet is encoded in 128 bits as a 2-element array of uint64_t - std::array encodedJet() { return encodedJet_; } - void setEncodedJet(std::array jet) { encodedJet_ = jet; } + // We store encodings both for Correlator internal usage and for Global Trigger + enum class HWEncoding { CT, GT }; + typedef std::array PackedJet; + const PackedJet& encodedJet(const HWEncoding encoding = HWEncoding::GT) const { + return encodedJet_[static_cast(encoding)]; + } + void setEncodedJet(const HWEncoding encoding, const PackedJet jet) { + encodedJet_[static_cast(encoding)] = jet; + } + + // Accessors to HW objects with ap_* types from encoded words + l1gt::Jet getHWJetGT() const { return l1gt::Jet::unpack(encodedJet(HWEncoding::GT)); } + l1ct::Jet getHWJetCT() const { return l1ct::Jet::unpack(encodedJet(HWEncoding::CT)); } private: float rawPt_; Constituents constituents_; - std::array encodedJet_ = {{0, 0}}; + std::array encodedJet_ = {{{{0, 0}}, {{0, 0}}}}; }; typedef std::vector PFJetCollection; 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 f56e7b36c2495..be802bc540b6e 100644 --- a/DataFormats/L1TParticleFlow/interface/datatypes.h +++ b/DataFormats/L1TParticleFlow/interface/datatypes.h @@ -39,6 +39,12 @@ 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; + typedef ap_ufixed<10, 1, AP_TRN, AP_SAT> srrtot_t; + 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; @@ -149,6 +155,11 @@ namespace l1ct { constexpr float Z0_LSB = 0.05; constexpr float DXY_LSB = 0.05; constexpr float PUPPIW_LSB = 1.0 / 256; + constexpr float MEANZ_OFFSET = 320.; + 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(); } inline float floatPt(pt2_t pt2) { return pt2.to_float(); } @@ -164,6 +175,9 @@ 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() / SRRTOT_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; } inline dpt_t makeDPt(int dpt) { return ap_fixed<18, 16>(dpt) >> 2; } @@ -194,6 +208,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 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)); }; 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/gt_datatypes.h b/DataFormats/L1TParticleFlow/interface/gt_datatypes.h index 3991831545009..64c2596d346c9 100644 --- a/DataFormats/L1TParticleFlow/interface/gt_datatypes.h +++ b/DataFormats/L1TParticleFlow/interface/gt_datatypes.h @@ -30,7 +30,7 @@ namespace l1gt { typedef ap_uint<1> valid_t; // E/gamma fields - typedef ap_fixed<11, 9> iso_t; + typedef ap_ufixed<11, 9> iso_t; typedef ap_uint<4> egquality_t; // tau fields diff --git a/DataFormats/L1TParticleFlow/interface/layer1_emulator.h b/DataFormats/L1TParticleFlow/interface/layer1_emulator.h index a257ed2efbda6..fcb415d426017 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; bool read(std::fstream &from); @@ -48,7 +48,6 @@ namespace l1ct { TkObj::clear(); src = nullptr; hwChi2 = 0; - hwStubs = 0; simPt = 0; simCaloEta = 0; simCaloPhi = 0; @@ -196,6 +195,7 @@ namespace l1ct { const l1t::PFTrack *srcTrack; // we use an index to the standalone object needed to retrieve a Ref when putting int sta_idx; + float idScore; bool read(std::fstream &from); bool write(std::fstream &to) const; void clear() { @@ -203,6 +203,7 @@ namespace l1ct { srcCluster = nullptr; srcTrack = nullptr; sta_idx = -1; + idScore = -999; clearIsoVars(); } @@ -334,7 +335,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..48af69b993de1 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,16 @@ 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_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; @@ -49,8 +61,12 @@ 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) { HadCaloObj ret; unsigned int start = 0; @@ -59,8 +75,13 @@ 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; } + + inline ap_uint pack_slim() const { return pack()(BITWIDTH_SLIM - 1, 0); } }; inline void clear(HadCaloObj &c) { c.clear(); } @@ -70,10 +91,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 +109,9 @@ namespace l1ct { hwEta = 0; hwPhi = 0; hwEmID = 0; + hwSrrTot = 0; + hwMeanZ = 0; + hwHoe = 0; } int intPt() const { return Scales::intPt(hwPt); } @@ -95,8 +122,14 @@ 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_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; - static const int BITWIDTH = pt_t::width + eta_t::width + phi_t::width + pt_t::width + emid_t::width; inline ap_uint pack() const { ap_uint ret; unsigned int start = 0; @@ -105,6 +138,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,8 +151,14 @@ 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; } + + inline ap_uint pack_slim() const { return pack()(BITWIDTH_SLIM - 1, 0); } }; inline void clear(EmCaloObj &c) { c.clear(); } @@ -130,6 +172,11 @@ namespace l1ct { z0_t hwZ0; dxy_t hwDxy; tkquality_t hwQuality; + stub_t hwStubs; + redChi2Bin_t hwRedChi2RZ; // 4 bits + redChi2Bin_t hwRedChi2RPhi; // 4 bits + //FIXME: 3 bits would be enough + redChi2Bin_t hwRedChi2Bend; // 4 bits enum TkQuality { PFLOOSE = 1, PFTIGHT = 2 }; bool isPFLoose() const { return hwQuality[0]; } @@ -140,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; + 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; } @@ -156,6 +204,10 @@ namespace l1ct { hwDxy = 0; hwCharge = false; hwQuality = 0; + hwStubs = 0; + hwRedChi2RZ = 0; + hwRedChi2RPhi = 0; + hwRedChi2Bend = 0; } int intPt() const { return Scales::intPt(hwPt); } @@ -174,8 +226,12 @@ namespace l1ct { float floatZ0() const { return Scales::floatZ0(hwZ0); } float floatDxy() const { return Scales::floatDxy(hwDxy); } - 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; + 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; @@ -188,6 +244,10 @@ 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, hwRedChi2RZ); + pack_into_bits(ret, start, hwRedChi2RPhi); + pack_into_bits(ret, start, hwRedChi2Bend); return ret; } inline static TkObj unpack(const ap_uint &src) { @@ -202,8 +262,13 @@ 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.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/interface/sums.h b/DataFormats/L1TParticleFlow/interface/sums.h index 195566f61afbe..b91aeec3b93d6 100644 --- a/DataFormats/L1TParticleFlow/interface/sums.h +++ b/DataFormats/L1TParticleFlow/interface/sums.h @@ -52,7 +52,7 @@ namespace l1ct { sum.valid = (hwPt != 0) || (hwSumPt != 0); sum.vector_pt = CTtoGT_pt(hwPt); sum.vector_phi = CTtoGT_phi(hwPhi); - sum.scalar_pt = CTtoGT_phi(hwSumPt); + sum.scalar_pt = CTtoGT_pt(hwSumPt); return sum; } }; 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..22ea4d1ecc7ef 100644 --- a/DataFormats/L1TParticleFlow/src/classes_def.xml +++ b/DataFormats/L1TParticleFlow/src/classes_def.xml @@ -1,6 +1,7 @@ - + + @@ -19,7 +20,8 @@ - + + @@ -33,7 +35,8 @@ - + + @@ -67,4 +70,3 @@ - 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/DataFormats/L1Trigger/interface/TkJetWord.h b/DataFormats/L1Trigger/interface/TkJetWord.h index a4f12dd57b627..8f45482c4e355 100644 --- a/DataFormats/L1Trigger/interface/TkJetWord.h +++ b/DataFormats/L1Trigger/interface/TkJetWord.h @@ -1,8 +1,10 @@ + #ifndef FIRMWARE_TkJetWord_h #define FIRMWARE_TkJetWord_h -// Class to store the 96-bit TkJet word for L1 Track Trigger. +// Class to store the 128-bit TkJet word for L1 Track Trigger. // Author: Emily MacDonald, updated by Benjamin Radburn-Smith (September 2022) +// 2nd update: George Karathanasis (Oct 2022) #include #include @@ -10,17 +12,15 @@ #include #include #include +#include "DataFormats/L1TrackTrigger/interface/TTTrack_TrackWord.h" namespace l1t { class TkJetWord { public: // ----------constants, enums and typedefs --------- - int INTPHI_PI = 720; - int INTPHI_TWOPI = 2 * INTPHI_PI; - float INTPT_LSB = 1 >> 5; - float ETAPHI_LSB = M_PI / (1 << 12); - float Z0_LSB = 0.05; + static constexpr double MAX_Z0 = 30.; + static constexpr double MAX_ETA = 8.; enum TkJetBitWidths { kPtSize = 16, @@ -30,33 +30,38 @@ namespace l1t { kZ0Size = 10, kNtSize = 5, kXtSize = 4, - kUnassignedSize = 8, - kTkJetWordSize = kPtSize + kGlbEtaSize + kGlbPhiSize + kZ0Size + kNtSize + kXtSize + kUnassignedSize, + kDispFlagSize = 1, + kUnassignedSize = 65, + kTkJetWordSize = + kPtSize + kGlbEtaSize + kGlbPhiSize + kZ0Size + kNtSize + kXtSize + kDispFlagSize + kUnassignedSize, }; enum TkJetBitLocations { kPtLSB = 0, kPtMSB = kPtLSB + TkJetBitWidths::kPtSize - 1, - kGlbEtaLSB = kPtMSB + 1, - kGlbEtaMSB = kGlbEtaLSB + TkJetBitWidths::kGlbEtaSize - 1, - kGlbPhiLSB = kGlbEtaMSB + 1, + kGlbPhiLSB = kPtMSB + 1, kGlbPhiMSB = kGlbPhiLSB + TkJetBitWidths::kGlbPhiSize - 1, - kZ0LSB = kGlbPhiMSB + 1, + kGlbEtaLSB = kGlbPhiMSB + 1, + kGlbEtaMSB = kGlbEtaLSB + TkJetBitWidths::kGlbEtaSize - 1, + kZ0LSB = kGlbEtaMSB + 1, kZ0MSB = kZ0LSB + TkJetBitWidths::kZ0Size - 1, kNtLSB = kZ0MSB + 1, kNtMSB = kNtLSB + TkJetBitWidths::kNtSize - 1, kXtLSB = kNtMSB + 1, kXtMSB = kXtLSB + TkJetBitWidths::kXtSize - 1, - kUnassignedLSB = kXtMSB + 1, + kDispFlagLSB = kXtMSB + 1, + kDispFlagMSB = kDispFlagLSB + TkJetBitWidths::kDispFlagSize - 1, + kUnassignedLSB = kDispFlagMSB + 1, kUnassignedMSB = kUnassignedLSB + TkJetBitWidths::kUnassignedSize - 1, }; typedef ap_ufixed pt_t; typedef ap_int glbeta_t; typedef ap_int glbphi_t; - typedef ap_int z0_t; // 40cm / 0.1 - typedef ap_uint nt_t; //number of tracks - typedef ap_uint nx_t; //number of tracks with xbit = 1 + typedef ap_int z0_t; // 40cm / 0.1 + typedef ap_uint nt_t; //number of tracks + typedef ap_uint nx_t; //number of tracks with xbit = 1 + typedef ap_uint dispflag_t; typedef ap_uint tkjetunassigned_t; // Unassigned bits typedef std::bitset tkjetword_bs_t; typedef ap_uint tkjetword_t; @@ -64,7 +69,14 @@ namespace l1t { public: // ----------Constructors -------------------------- TkJetWord() {} - TkJetWord(pt_t pt, glbeta_t eta, glbphi_t phi, z0_t z0, nt_t nt, nx_t nx, tkjetunassigned_t unassigned); + TkJetWord(pt_t pt, + glbeta_t eta, + glbphi_t phi, + z0_t z0, + nt_t nt, + nx_t nx, + dispflag_t dispflag, + tkjetunassigned_t unassigned); ~TkJetWord() {} @@ -109,6 +121,12 @@ namespace l1t { ret.V = tkJetWord()(TkJetBitLocations::kXtMSB, TkJetBitLocations::kXtLSB); return ret; } + dispflag_t dispFlagWord() const { + dispflag_t ret; + ret.V = tkJetWord()(TkJetBitLocations::kDispFlagMSB, TkJetBitLocations::kDispFlagLSB); + return ret; + } + tkjetunassigned_t unassignedWord() const { return tkJetWord()(TkJetBitLocations::kUnassignedMSB, TkJetBitLocations::kUnassignedLSB); } @@ -122,20 +140,37 @@ namespace l1t { unsigned int z0Bits() const { return z0Word().to_uint(); } unsigned int ntBits() const { return ntWord().to_uint(); } unsigned int xtBits() const { return xtWord().to_uint(); } + unsigned int dispFlagBits() const { return dispFlagWord().to_uint(); } unsigned int unassignedBits() const { return unassignedWord().to_uint(); } // These functions return the unpacked and converted values // These functions return real numbers converted from the digitized quantities by unpacking the 64-bit vertex word float pt() const { return ptWord().to_float(); } - float glbeta() const { return glbEtaWord().to_float() * ETAPHI_LSB; } - float glbphi() const { return glbPhiWord().to_float() * ETAPHI_LSB; } - float z0() const { return z0Word().to_float() * Z0_LSB; } + float glbeta() const { + return unpackSignedValue( + glbEtaWord(), TkJetBitWidths::kGlbEtaSize, (MAX_ETA) / (1 << TkJetBitWidths::kGlbEtaSize)); + } + float glbphi() const { + return unpackSignedValue( + glbPhiWord(), TkJetBitWidths::kGlbPhiSize, (2. * std::abs(M_PI)) / (1 << TkJetBitWidths::kGlbPhiSize)); + } + float z0() const { + return unpackSignedValue(z0Word(), TkJetBitWidths::kZ0Size, MAX_Z0 / (1 << TkJetBitWidths::kZ0Size)); + } int nt() const { return (ap_ufixed(ntWord())).to_int(); } int xt() const { return (ap_ufixed(xtWord())).to_int(); } + int dispflag() const { return (ap_ufixed(dispFlagWord())).to_int(); } unsigned int unassigned() const { return unassignedWord().to_uint(); } // ----------member functions (setters) ------------ - void setTkJetWord(pt_t pt, glbeta_t eta, glbphi_t phi, z0_t z0, nt_t nt, nx_t nx, tkjetunassigned_t unassigned); + void setTkJetWord(pt_t pt, + glbeta_t eta, + glbphi_t phi, + z0_t z0, + nt_t nt, + nx_t nx, + dispflag_t dispflag, + tkjetunassigned_t unassigned); private: // ----------private member functions -------------- diff --git a/DataFormats/L1Trigger/src/TkJetWord.cc b/DataFormats/L1Trigger/src/TkJetWord.cc index 2464c4a2db941..87897de7b07b0 100644 --- a/DataFormats/L1Trigger/src/TkJetWord.cc +++ b/DataFormats/L1Trigger/src/TkJetWord.cc @@ -4,26 +4,39 @@ #include "DataFormats/L1Trigger/interface/TkJetWord.h" namespace l1t { - TkJetWord::TkJetWord(pt_t pt, glbeta_t eta, glbphi_t phi, z0_t z0, nt_t nt, nx_t nx, tkjetunassigned_t unassigned) { - setTkJetWord(pt, eta, phi, z0, nt, nx, unassigned); + TkJetWord::TkJetWord(pt_t pt, + glbeta_t eta, + glbphi_t phi, + z0_t z0, + nt_t nt, + nx_t nx, + dispflag_t dispflag, + tkjetunassigned_t unassigned) { + setTkJetWord(pt, eta, phi, z0, nt, nx, dispflag, unassigned); } - void TkJetWord::setTkJetWord( - pt_t pt, glbeta_t eta, glbphi_t phi, z0_t z0, nt_t nt, nx_t nx, tkjetunassigned_t unassigned) { + void TkJetWord::setTkJetWord(pt_t pt, + glbeta_t eta, + glbphi_t phi, + z0_t z0, + nt_t nt, + nx_t nx, + dispflag_t dispflag, + tkjetunassigned_t unassigned) { // pack the TkJet word unsigned int offset = 0; for (unsigned int b = offset; b < (offset + TkJetBitWidths::kPtSize); b++) { tkJetWord_.set(b, pt[b - offset]); } offset += TkJetBitWidths::kPtSize; - for (unsigned int b = offset; b < (offset + TkJetBitWidths::kGlbEtaSize); b++) { - tkJetWord_.set(b, eta[b - offset]); - } - offset += TkJetBitWidths::kGlbEtaSize; for (unsigned int b = offset; b < (offset + TkJetBitWidths::kGlbPhiSize); b++) { tkJetWord_.set(b, phi[b - offset]); } offset += TkJetBitWidths::kGlbPhiSize; + for (unsigned int b = offset; b < (offset + TkJetBitWidths::kGlbEtaSize); b++) { + tkJetWord_.set(b, eta[b - offset]); + } + offset += TkJetBitWidths::kGlbEtaSize; for (unsigned int b = offset; b < (offset + TkJetBitWidths::kZ0Size); b++) { tkJetWord_.set(b, z0[b - offset]); } @@ -36,9 +49,13 @@ namespace l1t { tkJetWord_.set(b, nx[b - offset]); } offset += TkJetBitWidths::kXtSize; + for (unsigned int b = offset; b < (offset + TkJetBitWidths::kDispFlagSize); b++) { + tkJetWord_.set(b, nx[b - offset]); + } + offset += TkJetBitWidths::kDispFlagSize; for (unsigned int b = offset; b < (offset + TkJetBitWidths::kUnassignedSize); b++) { tkJetWord_.set(b, unassigned[b - offset]); } } -} //namespace l1t \ No newline at end of file +} //namespace l1t diff --git a/DataFormats/L1Trigger/src/classes_def.xml b/DataFormats/L1Trigger/src/classes_def.xml index 6b4f9444aaf04..b57276c8ebd03 100644 --- a/DataFormats/L1Trigger/src/classes_def.xml +++ b/DataFormats/L1Trigger/src/classes_def.xml @@ -304,7 +304,8 @@ - + + diff --git a/DataFormats/StdDictionaries/src/classes_def_others.xml b/DataFormats/StdDictionaries/src/classes_def_others.xml index e2c769fea5a1d..962d89e945c9b 100644 --- a/DataFormats/StdDictionaries/src/classes_def_others.xml +++ b/DataFormats/StdDictionaries/src/classes_def_others.xml @@ -12,6 +12,7 @@ + diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/l1t1PFPuppiJet70offMaxEta2p4_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/l1t1PFPuppiJet70offMaxEta2p4_cfi.py index 3eeda4877538e..af21009718979 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/modules/l1t1PFPuppiJet70offMaxEta2p4_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/l1t1PFPuppiJet70offMaxEta2p4_cfi.py @@ -10,5 +10,5 @@ endcap = cms.vdouble(42.4039, 1.33052, 0), overlap = cms.vdouble(24.8375, 1.4152, 0) ), - inputTag = cms.InputTag("l1tPhase1JetCalibrator9x9","Phase1L1TJetFromPfCandidates") + inputTag = cms.InputTag("l1tPhase1JetCalibrator9x9trimmed","Phase1L1TJetFromPfCandidates") ) diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/l1t2PFPuppiJet55offMaxEta2p4_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/l1t2PFPuppiJet55offMaxEta2p4_cfi.py index 9d17f392e98a7..4ce4aad79115b 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/modules/l1t2PFPuppiJet55offMaxEta2p4_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/l1t2PFPuppiJet55offMaxEta2p4_cfi.py @@ -10,5 +10,5 @@ endcap = cms.vdouble(42.4039, 1.33052, 0), overlap = cms.vdouble(24.8375, 1.4152, 0) ), - inputTag = cms.InputTag("l1tPhase1JetCalibrator9x9","Phase1L1TJetFromPfCandidates") + inputTag = cms.InputTag("l1tPhase1JetCalibrator9x9trimmed","Phase1L1TJetFromPfCandidates") ) diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/l1t4PFPuppiJet25OnlineMaxEta2p4_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/l1t4PFPuppiJet25OnlineMaxEta2p4_cfi.py index 77cf6d8138274..2c4f2b65f3f3f 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/modules/l1t4PFPuppiJet25OnlineMaxEta2p4_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/l1t4PFPuppiJet25OnlineMaxEta2p4_cfi.py @@ -5,5 +5,5 @@ MinEta = cms.double(-2.4), MinN = cms.int32(4), MinPt = cms.double(25.0), - inputTag = cms.InputTag("l1tPhase1JetCalibrator9x9","Phase1L1TJetFromPfCandidates") + inputTag = cms.InputTag("l1tPhase1JetCalibrator9x9trimmed","Phase1L1TJetFromPfCandidates") ) diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/l1t4PFPuppiJet40offMaxEta2p4_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/l1t4PFPuppiJet40offMaxEta2p4_cfi.py index 9e0c55b29708a..0ba10ab1e6e92 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/modules/l1t4PFPuppiJet40offMaxEta2p4_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/l1t4PFPuppiJet40offMaxEta2p4_cfi.py @@ -10,5 +10,5 @@ endcap = cms.vdouble(42.4039, 1.33052, 0), overlap = cms.vdouble(24.8375, 1.4152, 0) ), - inputTag = cms.InputTag("l1tPhase1JetCalibrator9x9","Phase1L1TJetFromPfCandidates") + inputTag = cms.InputTag("l1tPhase1JetCalibrator9x9trimmed","Phase1L1TJetFromPfCandidates") ) diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/l1tDoublePFPuppiJet112offMaxEta2p4_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/l1tDoublePFPuppiJet112offMaxEta2p4_cfi.py index 8bc28b6c44b02..22bef3d0647f0 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/modules/l1tDoublePFPuppiJet112offMaxEta2p4_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/l1tDoublePFPuppiJet112offMaxEta2p4_cfi.py @@ -10,6 +10,6 @@ endcap = cms.vdouble(42.4039, 1.33052, 0), overlap = cms.vdouble(24.8375, 1.4152, 0) ), - inputTag = cms.InputTag("l1tPhase1JetCalibrator9x9","Phase1L1TJetFromPfCandidates"), + inputTag = cms.InputTag("l1tPhase1JetCalibrator9x9trimmed","Phase1L1TJetFromPfCandidates"), saveTags = cms.bool(True) ) diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/l1tDoublePFPuppiJets112offMaxDeta1p6_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/l1tDoublePFPuppiJets112offMaxDeta1p6_cfi.py index f8c6ea593444d..79b52b0ae52a4 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/modules/l1tDoublePFPuppiJets112offMaxDeta1p6_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/l1tDoublePFPuppiJets112offMaxDeta1p6_cfi.py @@ -14,8 +14,8 @@ MinPt = cms.double(0.0), inputTag1 = cms.InputTag("l1tDoublePFPuppiJet112offMaxEta2p4"), inputTag2 = cms.InputTag("l1tDoublePFPuppiJet112offMaxEta2p4"), - originTag1 = cms.VInputTag(cms.InputTag("l1tPhase1JetCalibrator9x9","Phase1L1TJetFromPfCandidates")), - originTag2 = cms.VInputTag(cms.InputTag("l1tPhase1JetCalibrator9x9","Phase1L1TJetFromPfCandidates")), + originTag1 = cms.VInputTag(cms.InputTag("l1tPhase1JetCalibrator9x9trimmed","Phase1L1TJetFromPfCandidates")), + originTag2 = cms.VInputTag(cms.InputTag("l1tPhase1JetCalibrator9x9trimmed","Phase1L1TJetFromPfCandidates")), saveTags = cms.bool(True), triggerType1 = cms.int32(-116), triggerType2 = cms.int32(-116) diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/l1tPFPuppiHT400offMaxEta2p4_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/l1tPFPuppiHT400offMaxEta2p4_cfi.py index 648da465f075f..d387a3640748b 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/modules/l1tPFPuppiHT400offMaxEta2p4_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/l1tPFPuppiHT400offMaxEta2p4_cfi.py @@ -6,5 +6,5 @@ theScalings = cms.vdouble(50.0182, 1.0961, 0) ), TypeOfSum = cms.string('HT'), - inputTag = cms.InputTag("l1tPhase1JetSumsProducer","Sums") + inputTag = cms.InputTag("l1tPhase1JetSumsProducer9x9trimmed","Sums") ) diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/l1tPFPuppiHT450off_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/l1tPFPuppiHT450off_cfi.py index 090b3d17a160b..7d8cad4520804 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/modules/l1tPFPuppiHT450off_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/l1tPFPuppiHT450off_cfi.py @@ -6,5 +6,5 @@ theScalings = cms.vdouble(50.0182, 1.0961, 0) ), TypeOfSum = cms.string('HT'), - inputTag = cms.InputTag("l1tPhase1JetSumsProducer","Sums") + inputTag = cms.InputTag("l1tPhase1JetSumsProducer9x9trimmed","Sums") ) diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/l1tPFPuppiHT_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/l1tPFPuppiHT_cfi.py index 0f7a28e9ea4d9..7c37e7576eab5 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/modules/l1tPFPuppiHT_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/l1tPFPuppiHT_cfi.py @@ -1,7 +1,7 @@ import FWCore.ParameterSet.Config as cms l1tPFPuppiHT = cms.EDProducer("HLTHtMhtProducer", - jetsLabel = cms.InputTag("l1tPhase1JetCalibrator9x9","Phase1L1TJetFromPfCandidates"), + jetsLabel = cms.InputTag("l1tPhase1JetCalibrator9x9trimmed","Phase1L1TJetFromPfCandidates"), maxEtaJetHt = cms.double(2.4), minPtJetHt = cms.double(30.0) ) diff --git a/HLTrigger/Configuration/python/HLT_75e33/modules/l1tSinglePFPuppiJet230off_cfi.py b/HLTrigger/Configuration/python/HLT_75e33/modules/l1tSinglePFPuppiJet230off_cfi.py index 27463576e938a..25694aea4859c 100644 --- a/HLTrigger/Configuration/python/HLT_75e33/modules/l1tSinglePFPuppiJet230off_cfi.py +++ b/HLTrigger/Configuration/python/HLT_75e33/modules/l1tSinglePFPuppiJet230off_cfi.py @@ -10,5 +10,5 @@ endcap = cms.vdouble(42.4039, 1.33052, 0), overlap = cms.vdouble(24.8375, 1.4152, 0) ), - inputTag = cms.InputTag("l1tPhase1JetCalibrator9x9","Phase1L1TJetFromPfCandidates") + inputTag = cms.InputTag("l1tPhase1JetCalibrator9x9trimmed","Phase1L1TJetFromPfCandidates") ) diff --git a/L1Trigger/Configuration/python/L1Trigger_EventContent_cff.py b/L1Trigger/Configuration/python/L1Trigger_EventContent_cff.py index 559144d3dfaf0..6399f341f7984 100644 --- a/L1Trigger/Configuration/python/L1Trigger_EventContent_cff.py +++ b/L1Trigger/Configuration/python/L1Trigger_EventContent_cff.py @@ -199,6 +199,8 @@ def _appendPhase2Digis(obj): 'keep *_l1tPFTracksFromL1TracksHGCal_*_*', 'keep *_l1tSCPFL1PuppiCorrectedEmulator_*_*', 'keep *_l1tSCPFL1PuppiCorrectedEmulatorMHT_*_*', + 'keep *_l1tSCPFL1PuppiExtendedCorrectedEmulator_*_*', + 'keep *_l1tSCPFL1PuppiExtendedCorrectedEmulatorMHT_*_*', 'keep *_l1tPhase1JetProducer9x9_*_*', 'keep *_l1tPhase1JetCalibrator9x9_*_*', 'keep *_l1tPhase1JetSumsProducer9x9_*_*', @@ -210,7 +212,11 @@ def _appendPhase2Digis(obj): 'keep *_l1tLayer1HGCalNoTK_*_*', 'keep *_l1tLayer1HF_*_*', 'keep *_l1tLayer1_*_*', + 'keep *_l1tLayer1BarrelExtended_*_*', + 'keep *_l1tLayer1HGCalExtended_*_*', + 'keep *_l1tLayer1Extended_*_*', 'keep *_l1tLayer1EG_*_*', + 'keep *_l1tLayer2EG_*_*', 'keep *_l1tMETPFProducer_*_*', 'keep *_l1tNNTauProducer_*_*', 'keep *_l1tNNTauProducerPuppi_*_*', diff --git a/L1Trigger/Configuration/python/SimL1Emulator_cff.py b/L1Trigger/Configuration/python/SimL1Emulator_cff.py index 17600bf2e63a9..5dc7ba09569f6 100644 --- a/L1Trigger/Configuration/python/SimL1Emulator_cff.py +++ b/L1Trigger/Configuration/python/SimL1Emulator_cff.py @@ -151,9 +151,9 @@ from L1Trigger.L1TTrackMatch.l1tTrackerEtMiss_cfi import * from L1Trigger.L1TTrackMatch.l1tTrackerHTMiss_cfi import * # make the input tags consistent with the choice L1VertexFinder above -l1tTrackJets.L1PVertexCollection = ("l1tVertexFinder", "l1vertices") +l1tTrackJets.L1PVertexInputTag = ("l1tVertexFinderEmulator","l1verticesEmulation") l1tTrackFastJets.L1PrimaryVertexTag = ("l1tVertexFinder", "l1vertices") -l1tTrackJetsExtended.L1PVertexCollection = ("l1tVertexFinder", "l1vertices") +l1tTrackJetsExtended.L1PVertexInputTag = ("l1tVertexFinderEmulator","l1verticesEmulation") #l1tTrackerEtMiss.L1VertexInputTag = ("l1tVertexFinder", "l1vertices") #l1tTrackerEtMissExtended.L1VertexInputTag = ("l1tVertexFinder", "l1vertices") diff --git a/L1Trigger/DemonstratorTools/src/codecs_tkjets.cc b/L1Trigger/DemonstratorTools/src/codecs_tkjets.cc index 75de4a5c52e0b..069f23f5107da 100644 --- a/L1Trigger/DemonstratorTools/src/codecs_tkjets.cc +++ b/L1Trigger/DemonstratorTools/src/codecs_tkjets.cc @@ -34,13 +34,19 @@ namespace l1t::demo::codecs { //if (not x.test(TkJetWord::kValidLSB)) // break; - TkJetWord j(TkJetWord::pt_t(frames[f](TkJetWord::kPtMSB, TkJetWord::kPtLSB)), - TkJetWord::glbeta_t(frames[f](TkJetWord::kGlbEtaMSB, TkJetWord::kGlbEtaLSB)), - TkJetWord::glbphi_t(frames[f](TkJetWord::kGlbPhiMSB, TkJetWord::kGlbPhiLSB)), - TkJetWord::z0_t(frames[f](TkJetWord::kZ0MSB, TkJetWord::kZ0LSB)), - TkJetWord::nt_t(frames[f](TkJetWord::kNtMSB, TkJetWord::kNtLSB)), - TkJetWord::nx_t(frames[f](TkJetWord::kXtMSB, TkJetWord::kXtLSB)), - TkJetWord::tkjetunassigned_t(frames[f](TkJetWord::kUnassignedMSB, TkJetWord::kUnassignedLSB))); + TkJetWord j( + TkJetWord::pt_t(frames[f](TkJetWord::TkJetBitLocations::kPtMSB, TkJetWord::TkJetBitLocations::kPtLSB)), + TkJetWord::glbphi_t( + frames[f](TkJetWord::TkJetBitLocations::kGlbPhiMSB, TkJetWord::TkJetBitLocations::kGlbPhiLSB)), + TkJetWord::glbeta_t( + frames[f](TkJetWord::TkJetBitLocations::kGlbEtaMSB, TkJetWord::TkJetBitLocations::kGlbEtaLSB)), + TkJetWord::z0_t(frames[f](TkJetWord::TkJetBitLocations::kZ0MSB, TkJetWord::TkJetBitLocations::kZ0LSB)), + TkJetWord::nt_t(frames[f](TkJetWord::TkJetBitLocations::kNtMSB, TkJetWord::TkJetBitLocations::kNtLSB)), + TkJetWord::nx_t(frames[f](TkJetWord::TkJetBitLocations::kXtMSB, TkJetWord::TkJetBitLocations::kXtLSB)), + TkJetWord::dispflag_t( + frames[f](TkJetWord::TkJetBitLocations::kDispFlagMSB, TkJetWord::TkJetBitLocations::kDispFlagLSB)), + TkJetWord::tkjetunassigned_t( + frames[f](TkJetWord::TkJetBitLocations::kUnassignedMSB, TkJetWord::TkJetBitLocations::kUnassignedLSB))); tkJets.push_back(j); } diff --git a/L1Trigger/L1TTrackMatch/interface/L1Clustering.h b/L1Trigger/L1TTrackMatch/interface/L1Clustering.h deleted file mode 100644 index f97a11e8431ae..0000000000000 --- a/L1Trigger/L1TTrackMatch/interface/L1Clustering.h +++ /dev/null @@ -1,252 +0,0 @@ -#ifndef L1Trigger_L1TTrackMatch_L1Clustering_HH -#define L1Trigger_L1TTrackMatch_L1Clustering_HH -#include -#include -#include -#include - -//Each individual box in the eta and phi dimension. -// Also used to store final cluster data for each zbin. -struct EtaPhiBin { - float pTtot = 0; - int numtracks = 0; - int numttrks = 0; - int numtdtrks = 0; - int numttdtrks = 0; - bool used = false; - float phi; //average phi value (halfway b/t min and max) - float eta; //average eta value - std::vector trackidx; -}; - -//store important information for plots -struct MaxZBin { - int znum; //Numbered from 0 to nzbins (16, 32, or 64) in order - int nclust; //number of clusters in this bin - float zbincenter; - std::vector clusters; //list of all the clusters in this bin - float ht; //sum of all cluster pTs--only the zbin with the maximum ht is stored -}; - -inline std::vector L1_clustering(EtaPhiBin *phislice, int etaBins_, float etaStep_) { - std::vector clusters; - // Find eta-phibin with maxpT, make center of cluster, add neighbors if not already used - int nclust = 0; - - // get tracks in eta bins in increasing eta order - for (int etabin = 0; etabin < etaBins_; ++etabin) { - float my_pt = 0, previousbin_pt = 0; //, nextbin_pt=0, next2bin_pt=0; - float nextbin_pt = 0, nextbin2_pt = 0; - - // skip (already) used tracks - if (phislice[etabin].used) - continue; - my_pt = phislice[etabin].pTtot; - if (my_pt == 0) - continue; - //get previous bin pT - if (etabin > 0 && !phislice[etabin - 1].used) - previousbin_pt = phislice[etabin - 1].pTtot; - - // get next bins pt - if (etabin < etaBins_ - 1 && !phislice[etabin + 1].used) { - nextbin_pt = phislice[etabin + 1].pTtot; - if (etabin < etaBins_ - 2 && !phislice[etabin + 2].used) { - nextbin2_pt = phislice[etabin + 2].pTtot; - } - } - // check if pT of current cluster is higher than neighbors - if (my_pt < previousbin_pt || my_pt <= nextbin_pt) { - // if unused pT in the left neighbor, spit it out as a cluster - if (previousbin_pt > 0) { - clusters.push_back(phislice[etabin - 1]); - phislice[etabin - 1].used = true; - nclust++; - } - continue; //if it is not the local max pT skip - } - // here reach only unused local max clusters - clusters.push_back(phislice[etabin]); - phislice[etabin].used = true; //if current bin a cluster - if (previousbin_pt > 0) { - clusters[nclust].pTtot += previousbin_pt; - clusters[nclust].numtracks += phislice[etabin - 1].numtracks; - clusters[nclust].numtdtrks += phislice[etabin - 1].numtdtrks; - for (unsigned int itrk = 0; itrk < phislice[etabin - 1].trackidx.size(); itrk++) - clusters[nclust].trackidx.push_back(phislice[etabin - 1].trackidx[itrk]); - } - - if (my_pt >= nextbin2_pt && nextbin_pt > 0) { - clusters[nclust].pTtot += nextbin_pt; - clusters[nclust].numtracks += phislice[etabin + 1].numtracks; - clusters[nclust].numtdtrks += phislice[etabin + 1].numtdtrks; - for (unsigned int itrk = 0; itrk < phislice[etabin + 1].trackidx.size(); itrk++) - clusters[nclust].trackidx.push_back(phislice[etabin + 1].trackidx[itrk]); - phislice[etabin + 1].used = true; - } - - nclust++; - - } // for each etabin - - // Merge close-by clusters - for (int m = 0; m < nclust - 1; ++m) { - if (std::abs(clusters[m + 1].eta - clusters[m].eta) < 1.5 * etaStep_) { - if (clusters[m + 1].pTtot > clusters[m].pTtot) { - clusters[m].eta = clusters[m + 1].eta; - } - clusters[m].pTtot += clusters[m + 1].pTtot; - clusters[m].numtracks += clusters[m + 1].numtracks; // total ntrk - clusters[m].numtdtrks += clusters[m + 1].numtdtrks; // total ndisp - for (unsigned int itrk = 0; itrk < clusters[m + 1].trackidx.size(); itrk++) - clusters[m].trackidx.push_back(clusters[m + 1].trackidx[itrk]); - - // if remove the merged cluster - all the others must be closer to 0 - for (int m1 = m + 1; m1 < nclust - 1; ++m1) { - clusters[m1] = clusters[m1 + 1]; - //clusters.erase(clusters.begin()+m1); - } - // clusters[m1] = clusters[m1 + 1]; - clusters.erase(clusters.begin() + nclust); - nclust--; - m = -1; - } // end if clusters neighbor in eta - } // end for (m) loop - - return clusters; -} - -inline void Fill_L2Cluster(EtaPhiBin &bin, float pt, int ntrk, int ndtrk, std::vector trkidx) { - bin.pTtot += pt; - bin.numtracks += ntrk; - bin.numtdtrks += ndtrk; - for (unsigned int itrk = 0; itrk < trkidx.size(); itrk++) - bin.trackidx.push_back(trkidx[itrk]); -} - -inline float DPhi(float phi1, float phi2) { - float x = phi1 - phi2; - if (x >= M_PI) - x -= 2 * M_PI; - if (x < -1 * M_PI) - x += 2 * M_PI; - return x; -} - -inline std::vector L2_clustering(std::vector> &L1clusters, - int phiBins_, - float phiStep_, - float etaStep_) { - std::vector clusters; - for (int phibin = 0; phibin < phiBins_; ++phibin) { //Find eta-phibin with highest pT - if (L1clusters[phibin].empty()) - continue; - - // sort L1 clusters max -> min - sort(L1clusters[phibin].begin(), L1clusters[phibin].end(), [](struct EtaPhiBin &a, struct EtaPhiBin &b) { - return a.pTtot > b.pTtot; - }); - for (unsigned int imax = 0; imax < L1clusters[phibin].size(); ++imax) { - if (L1clusters[phibin][imax].used) - continue; - float pt_current = L1clusters[phibin][imax].pTtot; //current cluster (pt0) - float pt_next = 0; // next phi bin (pt1) - float pt_next2 = 0; // next to next phi bin2 (pt2) - int trk1 = 0; - int trk2 = 0; - int tdtrk1 = 0; - int tdtrk2 = 0; - std::vector trkidx1; - std::vector trkidx2; - clusters.push_back(L1clusters[phibin][imax]); - - L1clusters[phibin][imax].used = true; - - // if we are in the last phi bin, dont check phi+1 phi+2 - if (phibin == phiBins_ - 1) - continue; - std::vector used_already; //keep phi+1 clusters that have been used - for (unsigned int icluster = 0; icluster < L1clusters[phibin + 1].size(); ++icluster) { - if (L1clusters[phibin + 1][icluster].used) - continue; - if (std::abs(L1clusters[phibin + 1][icluster].eta - L1clusters[phibin][imax].eta) > 1.5 * etaStep_) - continue; - pt_next += L1clusters[phibin + 1][icluster].pTtot; - trk1 += L1clusters[phibin + 1][icluster].numtracks; - tdtrk1 += L1clusters[phibin + 1][icluster].numtdtrks; - for (unsigned int itrk = 0; itrk < L1clusters[phibin + 1][icluster].trackidx.size(); itrk++) - trkidx1.push_back(L1clusters[phibin + 1][icluster].trackidx[itrk]); - used_already.push_back(icluster); - } - - if (pt_next < pt_current) { // if pt1 used_already2; //keep used clusters in phi+2 - for (unsigned int icluster = 0; icluster < L1clusters[phibin + 2].size(); ++icluster) { - if (L1clusters[phibin + 2][icluster].used) - continue; - if (std::abs(L1clusters[phibin + 2][icluster].eta - L1clusters[phibin][imax].eta) > 1.5 * etaStep_) - continue; - pt_next2 += L1clusters[phibin + 2][icluster].pTtot; - trk2 += L1clusters[phibin + 2][icluster].numtracks; - tdtrk2 += L1clusters[phibin + 2][icluster].numtdtrks; - for (unsigned int itrk = 0; itrk < L1clusters[phibin + 2][icluster].trackidx.size(); itrk++) - trkidx2.push_back(L1clusters[phibin + 2][icluster].trackidx[itrk]); - used_already2.push_back(icluster); - } - if (pt_next2 < pt_next) { - std::vector trkidx_both; - trkidx_both.reserve(trkidx1.size() + trkidx2.size()); - trkidx_both.insert(trkidx_both.end(), trkidx1.begin(), trkidx1.end()); - trkidx_both.insert(trkidx_both.end(), trkidx2.begin(), trkidx2.end()); - Fill_L2Cluster(clusters[clusters.size() - 1], pt_next + pt_next2, trk1 + trk2, tdtrk1 + tdtrk2, trkidx_both); - clusters[clusters.size() - 1].phi = L1clusters[phibin + 1][used_already[0]].phi; - for (unsigned int iused : used_already) - L1clusters[phibin + 1][iused].used = true; - for (unsigned int iused : used_already2) - L1clusters[phibin + 2][iused].used = true; - } - } // for each L1 cluster - } // for each phibin - - int nclust = clusters.size(); - - // merge close-by clusters - for (int m = 0; m < nclust - 1; ++m) { - for (int n = m + 1; n < nclust; ++n) { - if (clusters[n].eta != clusters[m].eta) - continue; - if (std::abs(DPhi(clusters[n].phi, clusters[m].phi)) > 1.5 * phiStep_) - continue; - - if (clusters[n].pTtot > clusters[m].pTtot) - clusters[m].phi = clusters[n].phi; - - clusters[m].pTtot += clusters[n].pTtot; - clusters[m].numtracks += clusters[n].numtracks; - clusters[m].numtdtrks += clusters[n].numtdtrks; - for (unsigned int itrk = 0; itrk < clusters[n].trackidx.size(); itrk++) - clusters[m].trackidx.push_back(clusters[n].trackidx[itrk]); - for (int m1 = n; m1 < nclust - 1; ++m1) - clusters[m1] = clusters[m1 + 1]; - clusters.erase(clusters.begin() + nclust); - - nclust--; - m = -1; - } // end of n-loop - } // end of m-loop - return clusters; -} -#endif diff --git a/L1Trigger/L1TTrackMatch/interface/L1TrackJetEmulationProducer.h b/L1Trigger/L1TTrackMatch/interface/L1TrackJetEmulationProducer.h deleted file mode 100644 index dd46370d13af4..0000000000000 --- a/L1Trigger/L1TTrackMatch/interface/L1TrackJetEmulationProducer.h +++ /dev/null @@ -1,79 +0,0 @@ -#ifndef L1Trigger_L1TTrackMatch_L1TrackJetEmulationProducer_HH -#define L1Trigger_L1TTrackMatch_L1TrackJetEmulationProducer_HH -#include -#include -#include -#include -#include -#include -#include "DataFormats/L1Trigger/interface/TkJetWord.h" - -//For precision studies -const int PT_EXTRABITS = 0; -const int ETA_EXTRABITS = 0; -const int PHI_EXTRABITS = 0; -const int Z0_EXTRABITS = 0; - -typedef ap_ufixed<16 + PT_EXTRABITS, 11, AP_TRN, AP_SAT> pt_intern; -typedef ap_int<14 + ETA_EXTRABITS> glbeta_intern; -typedef ap_int<14 + PHI_EXTRABITS> glbphi_intern; -typedef ap_int<10 + Z0_EXTRABITS> z0_intern; // 40cm / 0.1 - -namespace convert { - const int INTPHI_PI = 720; - const int INTPHI_TWOPI = 2 * INTPHI_PI; - static const float INTPT_LSB_POW = pow(2.0, -5 - PT_EXTRABITS); - static const float INTPT_LSB = INTPT_LSB_POW; - static const float ETA_LSB_POW = pow(2.0, -1 * ETA_EXTRABITS); - static const float ETA_LSB = M_PI / pow(2.0, 12) * ETA_LSB_POW; - static const float PHI_LSB_POW = pow(2.0, -1 * PHI_EXTRABITS); - static const float PHI_LSB = M_PI / pow(2.0, 12) * PHI_LSB_POW; - static const float Z0_LSB_POW = pow(2.0, -1 * Z0_EXTRABITS); - static const float Z0_LSB = 0.05 * Z0_LSB_POW; - inline float floatPt(pt_intern pt) { return pt.to_float(); } - inline int intPt(pt_intern pt) { return (ap_ufixed<18 + PT_EXTRABITS, 13 + PT_EXTRABITS>(pt)).to_int(); } - inline float floatEta(glbeta_intern eta) { return eta.to_float() * ETA_LSB; } - inline float floatPhi(glbphi_intern phi) { return phi.to_float() * PHI_LSB; } - inline float floatZ0(z0_intern z0) { return z0.to_float() * Z0_LSB; } - - inline pt_intern makePt(int pt) { return ap_ufixed<18 + PT_EXTRABITS, 13 + PT_EXTRABITS>(pt); } - inline pt_intern makePtFromFloat(float pt) { return pt_intern(INTPT_LSB_POW * round(pt / INTPT_LSB_POW)); } - inline z0_intern makeZ0(float z0) { return z0_intern(round(z0 / Z0_LSB)); } - - inline ap_uint ptToInt(pt_intern pt) { - // note: this can be synthethized, e.g. when pT is used as intex in a LUT - ap_uint ret = 0; - ret(pt_intern::width - 1, 0) = pt(pt_intern::width - 1, 0); - return ret; - } - - inline glbeta_intern makeGlbEta(float eta) { return round(eta / ETA_LSB); } - inline glbeta_intern makeGlbEtaRoundEven(float eta) { - glbeta_intern ghweta = round(eta / ETA_LSB); - return (ghweta % 2) ? glbeta_intern(ghweta + 1) : ghweta; - } - - inline glbphi_intern makeGlbPhi(float phi) { return round(phi / PHI_LSB); } - -}; // namespace convert - -//Each individual box in the eta and phi dimension. -// Also used to store final cluster data for each zbin. -struct TrackJetEmulationEtaPhiBin { - pt_intern pTtot; - l1t::TkJetWord::nt_t ntracks; - l1t::TkJetWord::nx_t nxtracks; - bool used; - glbphi_intern phi; //average phi value (halfway b/t min and max) - glbeta_intern eta; //average eta value -}; - -//store important information for plots -struct TrackJetEmulationMaxZBin { - int znum; //Numbered from 0 to nzbins (16, 32, or 64) in order - int nclust; //number of clusters in this bin - z0_intern zbincenter; - TrackJetEmulationEtaPhiBin *clusters; //list of all the clusters in this bin - pt_intern ht; //sum of all cluster pTs--only the zbin with the maximum ht is stored -}; -#endif diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetClustering.h b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetClustering.h new file mode 100644 index 0000000000000..c4fdf3dfd1c02 --- /dev/null +++ b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetClustering.h @@ -0,0 +1,365 @@ +#ifndef L1Trigger_L1TTrackMatch_L1TrackJetClustering_HH +#define L1Trigger_L1TTrackMatch_L1TrackJetClustering_HH +#include +#include +#include +#include +#include +#include +#include "DataFormats/L1Trigger/interface/TkJetWord.h" +#include "DataFormats/L1TrackTrigger/interface/TTTrack_TrackWord.h" +#include "DataFormats/L1TrackTrigger/interface/TTTypes.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +namespace l1ttrackjet { + //For precision studies + const unsigned int PT_INTPART_BITS{9}; + const unsigned int ETA_INTPART_BITS{3}; + const unsigned int kExtraGlobalPhiBit{4}; + + typedef ap_ufixed pt_intern; + typedef ap_fixed glbeta_intern; + typedef ap_int glbphi_intern; + typedef ap_int z0_intern; // 40cm / 0.1 + typedef ap_uint d0_intern; + + inline const unsigned int DoubleToBit(double value, unsigned int maxBits, double step) { + unsigned int digitized_value = std::floor(std::abs(value) / step); + unsigned int digitized_maximum = (1 << (maxBits - 1)) - 1; // The remove 1 bit from nBits to account for the sign + if (digitized_value > digitized_maximum) + digitized_value = digitized_maximum; + if (value < 0) + digitized_value = (1 << maxBits) - digitized_value; // two's complement encoding + return digitized_value; + } + inline const double BitToDouble(unsigned int bits, unsigned int maxBits, double step) { + int isign = 1; + unsigned int digitized_maximum = (1 << maxBits) - 1; + if (bits & (1 << (maxBits - 1))) { // check the sign + isign = -1; + bits = (1 << (maxBits + 1)) - bits; + } + return (double(bits & digitized_maximum) + 0.5) * step * isign; + } + + // eta/phi clusters - simulation + struct EtaPhiBin { + float pTtot; + int ntracks; + int nxtracks; + bool used; + float phi; //average phi value (halfway b/t min and max) + float eta; //average eta value + std::vector trackidx; + }; + // z bin struct - simulation (used if z bin are many) + struct MaxZBin { + int znum; //Numbered from 0 to nzbins (16, 32, or 64) in order + int nclust; //number of clusters in this bin + float zbincenter; + std::vector clusters; //list of all the clusters in this bin + float ht; //sum of all cluster pTs--only the zbin with the maximum ht is stored + }; + + // eta/phi clusters - emulation + struct TrackJetEmulationEtaPhiBin { + pt_intern pTtot; + l1t::TkJetWord::nt_t ntracks; + l1t::TkJetWord::nx_t nxtracks; + bool used; + glbphi_intern phi; //average phi value (halfway b/t min and max) + glbeta_intern eta; //average eta value + std::vector trackidx; + }; + + // z bin struct - emulation (used if z bin are many) + struct TrackJetEmulationMaxZBin { + int znum; //Numbered from 0 to nzbins (16, 32, or 64) in order + int nclust; //number of clusters in this bin + z0_intern zbincenter; + std::vector clusters; //list of all the clusters in this bin + pt_intern ht; //sum of all cluster pTs--only the zbin with the maximum ht is stored + }; + + // track quality cuts + inline bool TrackQualitySelection(int trk_nstub, + double trk_chi2, + double trk_bendchi2, + double nStubs4PromptBend_, + double nStubs5PromptBend_, + double nStubs4PromptChi2_, + double nStubs5PromptChi2_, + double nStubs4DisplacedBend_, + double nStubs5DisplacedBend_, + double nStubs4DisplacedChi2_, + double nStubs5DisplacedChi2_, + bool displaced_) { + bool PassQuality = false; + if (!displaced_) { + if (trk_nstub == 4 && trk_bendchi2 < nStubs4PromptBend_ && + trk_chi2 < nStubs4PromptChi2_) // 4 stubs are the lowest track quality and have different cuts + PassQuality = true; + if (trk_nstub > 4 && trk_bendchi2 < nStubs5PromptBend_ && + trk_chi2 < nStubs5PromptChi2_) // above 4 stubs diffent selection imposed (genrally looser) + PassQuality = true; + } else { + if (trk_nstub == 4 && trk_bendchi2 < nStubs4DisplacedBend_ && + trk_chi2 < nStubs4DisplacedChi2_) // 4 stubs are the lowest track quality and have different cuts + PassQuality = true; + if (trk_nstub > 4 && trk_bendchi2 < nStubs5DisplacedBend_ && + trk_chi2 < nStubs5DisplacedChi2_) // above 4 stubs diffent selection imposed (genrally looser) + PassQuality = true; + } + return PassQuality; + } + + // L1 clustering (in eta) + template + inline std::vector L1_clustering(T *phislice, int etaBins_, Eta etaStep_) { + std::vector clusters; + // Find eta bin with maxpT, make center of cluster, add neighbors if not already used + int nclust = 0; + + // get tracks in eta bins in increasing eta order + for (int etabin = 0; etabin < etaBins_; ++etabin) { + Pt my_pt = 0; + Pt previousbin_pt = 0; + Pt nextbin_pt = 0; + Pt nextbin2_pt = 0; + + // skip (already) used tracks + if (phislice[etabin].used) + continue; + + my_pt = phislice[etabin].pTtot; + if (my_pt == 0) + continue; + + //get previous bin pT + if (etabin > 0 && !phislice[etabin - 1].used) + previousbin_pt = phislice[etabin - 1].pTtot; + + // get next bins pt + if (etabin < etaBins_ - 1 && !phislice[etabin + 1].used) { + nextbin_pt = phislice[etabin + 1].pTtot; + if (etabin < etaBins_ - 2 && !phislice[etabin + 2].used) { + nextbin2_pt = phislice[etabin + 2].pTtot; + } + } + // check if pT of current cluster is higher than neighbors + if (my_pt < previousbin_pt || my_pt <= nextbin_pt) { + // if unused pT in the left neighbor, spit it out as a cluster + if (previousbin_pt > 0) { + clusters.push_back(phislice[etabin - 1]); + phislice[etabin - 1].used = true; + nclust++; + } + continue; //if it is not the local max pT skip + } + // here reach only unused local max clusters + clusters.push_back(phislice[etabin]); + phislice[etabin].used = true; //if current bin a cluster + if (previousbin_pt > 0) { + clusters[nclust].pTtot += previousbin_pt; + clusters[nclust].ntracks += phislice[etabin - 1].ntracks; + clusters[nclust].nxtracks += phislice[etabin - 1].nxtracks; + for (unsigned int itrk = 0; itrk < phislice[etabin - 1].trackidx.size(); itrk++) + clusters[nclust].trackidx.push_back(phislice[etabin - 1].trackidx[itrk]); + } + + if (my_pt >= nextbin2_pt && nextbin_pt > 0) { + clusters[nclust].pTtot += nextbin_pt; + clusters[nclust].ntracks += phislice[etabin + 1].ntracks; + clusters[nclust].nxtracks += phislice[etabin + 1].nxtracks; + for (unsigned int itrk = 0; itrk < phislice[etabin + 1].trackidx.size(); itrk++) + clusters[nclust].trackidx.push_back(phislice[etabin + 1].trackidx[itrk]); + phislice[etabin + 1].used = true; + } + + nclust++; + + } // for each etabin + + // Merge close-by clusters + for (int m = 0; m < nclust - 1; ++m) { + if (((clusters[m + 1].eta - clusters[m].eta) < (3 * etaStep_) / 2) && + (-(3 * etaStep_) / 2 < (clusters[m + 1].eta - clusters[m].eta))) { + if (clusters[m + 1].pTtot > clusters[m].pTtot) { + clusters[m].eta = clusters[m + 1].eta; + } + clusters[m].pTtot += clusters[m + 1].pTtot; + clusters[m].ntracks += clusters[m + 1].ntracks; // total ntrk + clusters[m].nxtracks += clusters[m + 1].nxtracks; // total ndisp + for (unsigned int itrk = 0; itrk < clusters[m + 1].trackidx.size(); itrk++) + clusters[m].trackidx.push_back(clusters[m + 1].trackidx[itrk]); + + // if remove the merged cluster - all the others must be closer to 0 + for (int m1 = m + 1; m1 < nclust - 1; ++m1) + clusters[m1] = clusters[m1 + 1]; + + clusters.erase(clusters.begin() + nclust); + nclust--; + } // end if for cluster merging + } // end for (m) loop + + return clusters; + } + + // Fill L2 clusters (helper function) + template + inline void Fill_L2Cluster(T &bin, Pt pt, int ntrk, int ndtrk, std::vector trkidx) { + bin.pTtot += pt; + bin.ntracks += ntrk; + bin.nxtracks += ndtrk; + for (unsigned int itrk = 0; itrk < trkidx.size(); itrk++) + bin.trackidx.push_back(trkidx[itrk]); + } + + inline glbphi_intern DPhi(glbphi_intern phi1, glbphi_intern phi2) { + glbphi_intern x = phi1 - phi2; + if (x >= DoubleToBit( + M_PI, TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, TTTrack_TrackWord::stepPhi0)) + x -= DoubleToBit( + 2 * M_PI, TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, TTTrack_TrackWord::stepPhi0); + if (x < DoubleToBit(-1 * M_PI, + TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, + TTTrack_TrackWord::stepPhi0)) + x += DoubleToBit( + 2 * M_PI, TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, TTTrack_TrackWord::stepPhi0); + return x; + } + + inline float DPhi(float phi1, float phi2) { + float x = phi1 - phi2; + if (x >= M_PI) + x -= 2 * M_PI; + if (x < -1 * M_PI) + x += 2 * M_PI; + return x; + } + + // L2 clustering (in phi) + template + inline std::vector L2_clustering(std::vector > &L1clusters, + int phiBins_, + Phi phiStep_, + Eta etaStep_) { + std::vector clusters; + for (int phibin = 0; phibin < phiBins_; ++phibin) { //Find eta-phibin with highest pT + if (L1clusters[phibin].empty()) + continue; + + // sort L1 clusters max -> min + sort(L1clusters[phibin].begin(), L1clusters[phibin].end(), [](T &a, T &b) { return a.pTtot > b.pTtot; }); + + for (unsigned int imax = 0; imax < L1clusters[phibin].size(); ++imax) { + if (L1clusters[phibin][imax].used) + continue; + Pt pt_current = L1clusters[phibin][imax].pTtot; //current cluster (pt0) + Pt pt_next = 0; // next phi bin (pt1) + Pt pt_next2 = 0; // next to next phi bin2 (pt2) + int trk1 = 0; + int trk2 = 0; + int tdtrk1 = 0; + int tdtrk2 = 0; + std::vector trkidx1; + std::vector trkidx2; + clusters.push_back(L1clusters[phibin][imax]); + + L1clusters[phibin][imax].used = true; + + // if we are in the last phi bin, dont check phi+1 phi+2 + if (phibin == phiBins_ - 1) + continue; + + std::vector used_already; //keep phi+1 clusters that have been used + for (unsigned int icluster = 0; icluster < L1clusters[phibin + 1].size(); ++icluster) { + if (L1clusters[phibin + 1][icluster].used) + continue; + + if (((L1clusters[phibin + 1][icluster].eta - L1clusters[phibin][imax].eta) > (3 * etaStep_) / 2) || + ((L1clusters[phibin + 1][icluster].eta - L1clusters[phibin][imax].eta) < -(3 * etaStep_) / 2)) + continue; + + pt_next += L1clusters[phibin + 1][icluster].pTtot; + trk1 += L1clusters[phibin + 1][icluster].ntracks; + tdtrk1 += L1clusters[phibin + 1][icluster].nxtracks; + + for (unsigned int itrk = 0; itrk < L1clusters[phibin + 1][icluster].trackidx.size(); itrk++) + trkidx1.push_back(L1clusters[phibin + 1][icluster].trackidx[itrk]); + used_already.push_back(icluster); + } + + if (pt_next < pt_current) { // if pt1(clusters[clusters.size() - 1], pt_next, trk1, tdtrk1, trkidx1); + for (unsigned int iused : used_already) + L1clusters[phibin + 1][iused].used = true; + continue; + } + // if phi = next to last bin there is no "next to next" + if (phibin == phiBins_ - 2) { + Fill_L2Cluster(clusters[clusters.size() - 1], pt_next, trk1, tdtrk1, trkidx1); + clusters[clusters.size() - 1].phi = L1clusters[phibin + 1][used_already[0]].phi; + for (unsigned int iused : used_already) + L1clusters[phibin + 1][iused].used = true; + continue; + } + std::vector used_already2; //keep used clusters in phi+2 + for (unsigned int icluster = 0; icluster < L1clusters[phibin + 2].size(); ++icluster) { + if (L1clusters[phibin + 2][icluster].used) + continue; + if (((L1clusters[phibin + 2][icluster].eta - L1clusters[phibin][imax].eta) > (3 * etaStep_) / 2) || + ((L1clusters[phibin + 2][icluster].eta - L1clusters[phibin][imax].eta) < -(3 * etaStep_) / 2)) + continue; + pt_next2 += L1clusters[phibin + 2][icluster].pTtot; + trk2 += L1clusters[phibin + 2][icluster].ntracks; + tdtrk2 += L1clusters[phibin + 2][icluster].nxtracks; + + for (unsigned int itrk = 0; itrk < L1clusters[phibin + 2][icluster].trackidx.size(); itrk++) + trkidx2.push_back(L1clusters[phibin + 2][icluster].trackidx[itrk]); + used_already2.push_back(icluster); + } + if (pt_next2 < pt_next) { + std::vector trkidx_both; + trkidx_both.reserve(trkidx1.size() + trkidx2.size()); + trkidx_both.insert(trkidx_both.end(), trkidx1.begin(), trkidx1.end()); + trkidx_both.insert(trkidx_both.end(), trkidx2.begin(), trkidx2.end()); + Fill_L2Cluster( + clusters[clusters.size() - 1], pt_next + pt_next2, trk1 + trk2, tdtrk1 + tdtrk2, trkidx_both); + clusters[clusters.size() - 1].phi = L1clusters[phibin + 1][used_already[0]].phi; + for (unsigned int iused : used_already) + L1clusters[phibin + 1][iused].used = true; + for (unsigned int iused : used_already2) + L1clusters[phibin + 2][iused].used = true; + } + } // for each L1 cluster + } // for each phibin + + int nclust = clusters.size(); + + // merge close-by clusters + for (int m = 0; m < nclust - 1; ++m) { + for (int n = m + 1; n < nclust; ++n) { + if (clusters[n].eta != clusters[m].eta) + continue; + if ((DPhi(clusters[n].phi, clusters[m].phi) > (3 * phiStep_) / 2) || + (DPhi(clusters[n].phi, clusters[m].phi) < -(3 * phiStep_) / 2)) + continue; + if (clusters[n].pTtot > clusters[m].pTtot) + clusters[m].phi = clusters[n].phi; + clusters[m].pTtot += clusters[n].pTtot; + clusters[m].ntracks += clusters[n].ntracks; + clusters[m].nxtracks += clusters[n].nxtracks; + for (unsigned int itrk = 0; itrk < clusters[n].trackidx.size(); itrk++) + clusters[m].trackidx.push_back(clusters[n].trackidx[itrk]); + for (int m1 = n; m1 < nclust - 1; ++m1) + clusters[m1] = clusters[m1 + 1]; + clusters.erase(clusters.begin() + nclust); + + nclust--; + } // end of n-loop + } // end of m-loop + return clusters; + } +} // namespace l1ttrackjet +#endif diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetEmulationProducer.cc b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetEmulationProducer.cc deleted file mode 100644 index 701be4b196efa..0000000000000 --- a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetEmulationProducer.cc +++ /dev/null @@ -1,681 +0,0 @@ -/* Software to emulate the hardware 2-layer jet-finding algorithm (fixed-point). *Layers 1 and 2* - * - * 2021 - * - * Created based on Rishi Patel's L1TrackJetProducer.cc file. - * Authors: Samuel Edwin Leigh, Tyler Wu - * Rutgers, the State University of New Jersey - * Revolutionary for 250 years - */ - -//Holds data from tracks, converted from their integer versions. - -// system include files - -#include "DataFormats/Common/interface/Ref.h" -#include "DataFormats/L1TCorrelator/interface/TkJet.h" -#include "DataFormats/L1TCorrelator/interface/TkJetFwd.h" -#include "DataFormats/L1TrackTrigger/interface/TTTypes.h" -#include "DataFormats/L1TrackTrigger/interface/TTTrack.h" -#include "DataFormats/L1TrackTrigger/interface/TTTrack_TrackWord.h" -#include "DataFormats/L1Trigger/interface/TkJetWord.h" -#include "DataFormats/L1Trigger/interface/VertexWord.h" -// user include files -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/stream/EDProducer.h" -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/EventSetup.h" -#include "FWCore/Framework/interface/MakerMacros.h" -#include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "FWCore/Utilities/interface/StreamID.h" -#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" -#include "Geometry/CommonTopologies/interface/PixelGeomDetUnit.h" -#include "Geometry/CommonTopologies/interface/PixelGeomDetType.h" -#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" - -#include "DataFormats/L1Trigger/interface/Vertex.h" - -#include -#include -#include -#include -#include -#include -#include -#include "TH1D.h" -#include "TH2D.h" -#include -#include -#include "L1Trigger/L1TTrackMatch/interface/L1TrackJetEmulationProducer.h" - -using namespace std; -using namespace edm; -class L1TrackJetEmulationProducer : public stream::EDProducer<> { -public: - explicit L1TrackJetEmulationProducer(const ParameterSet &); - ~L1TrackJetEmulationProducer() override; - typedef TTTrack L1TTTrackType; - typedef vector L1TTTrackCollectionType; - - static void fillDescriptions(ConfigurationDescriptions &descriptions); - bool trackQualityCuts(int trk_nstub, float trk_chi2, float trk_bendchi2); - void L2_cluster(vector> L1TrkPtrs_, TrackJetEmulationMaxZBin &mzb); - virtual TrackJetEmulationEtaPhiBin *L1_cluster(TrackJetEmulationEtaPhiBin *phislice); - -private: - void beginStream(StreamID) override; - void produce(Event &, const EventSetup &) override; - void endStream() override; - - // ----------member data --------------------------- - - vector> L1TrkPtrs_; - vector zBinCount_; - vector tdtrk_; - const float MaxDzTrackPV; - const float trkZMax_; - const float trkPtMax_; - const float trkPtMin_; - const float trkEtaMax_; - const float trkChi2dofMax_; - const float trkBendChi2Max_; - const int trkNPSStubMin_; - const double minTrkJetpT_; - int etaBins_; - int phiBins_; - int zBins_; - float d0CutNStubs4_; - float d0CutNStubs5_; - int lowpTJetMinTrackMultiplicity_; - float lowpTJetMinpT_; - int highpTJetMinTrackMultiplicity_; - float highpTJetMinpT_; - bool displaced_; - float nStubs4DisplacedChi2_; - float nStubs5DisplacedChi2_; - float nStubs4DisplacedBend_; - float nStubs5DisplacedBend_; - int nDisplacedTracks_; - - float PVz; - z0_intern zStep_; - glbeta_intern etaStep_; - glbphi_intern phiStep_; - - const EDGetTokenT>> trackToken_; - const EDGetTokenT PVtxToken_; - ESGetToken tTopoToken_; -}; - -L1TrackJetEmulationProducer::L1TrackJetEmulationProducer(const ParameterSet &iConfig) - : MaxDzTrackPV((float)iConfig.getParameter("MaxDzTrackPV")), - trkZMax_((float)iConfig.getParameter("trk_zMax")), - trkPtMax_((float)iConfig.getParameter("trk_ptMax")), - trkPtMin_((float)iConfig.getParameter("trk_ptMin")), - trkEtaMax_((float)iConfig.getParameter("trk_etaMax")), - trkChi2dofMax_((float)iConfig.getParameter("trk_chi2dofMax")), - trkBendChi2Max_((float)iConfig.getParameter("trk_bendChi2Max")), - trkNPSStubMin_((int)iConfig.getParameter("trk_nPSStubMin")), - minTrkJetpT_(iConfig.getParameter("minTrkJetpT")), - etaBins_((int)iConfig.getParameter("etaBins")), - phiBins_((int)iConfig.getParameter("phiBins")), - zBins_((int)iConfig.getParameter("zBins")), - d0CutNStubs4_((float)iConfig.getParameter("d0_cutNStubs4")), - d0CutNStubs5_((float)iConfig.getParameter("d0_cutNStubs5")), - lowpTJetMinTrackMultiplicity_((int)iConfig.getParameter("lowpTJetMinTrackMultiplicity")), - lowpTJetMinpT_((float)iConfig.getParameter("lowpTJetMinpT")), - highpTJetMinTrackMultiplicity_((int)iConfig.getParameter("highpTJetMinTrackMultiplicity")), - highpTJetMinpT_((float)iConfig.getParameter("highpTJetMinpT")), - displaced_(iConfig.getParameter("displaced")), - nStubs4DisplacedChi2_((float)iConfig.getParameter("nStubs4DisplacedChi2")), - nStubs5DisplacedChi2_((float)iConfig.getParameter("nStubs5DisplacedChi2")), - nStubs4DisplacedBend_((float)iConfig.getParameter("nStubs4Displacedbend")), - nStubs5DisplacedBend_((float)iConfig.getParameter("nStubs5Displacedbend")), - nDisplacedTracks_((int)iConfig.getParameter("nDisplacedTracks")), - trackToken_(consumes>>(iConfig.getParameter("L1TrackInputTag"))), - PVtxToken_(consumes(iConfig.getParameter("VertexInputTag"))), - tTopoToken_(esConsumes(edm::ESInputTag("", ""))) { - zStep_ = convert::makeZ0(2.0 * trkZMax_ / (zBins_ + 1)); - etaStep_ = convert::makeGlbEta(2.0 * trkEtaMax_ / etaBins_); //etaStep is the width of an etabin - phiStep_ = convert::makeGlbPhi(2.0 * M_PI / phiBins_); ////phiStep is the width of a phibin - - if (displaced_) - produces("L1TrackJetsExtended"); - else - produces("L1TrackJets"); -} - -L1TrackJetEmulationProducer::~L1TrackJetEmulationProducer() {} - -void L1TrackJetEmulationProducer::produce(Event &iEvent, const EventSetup &iSetup) { - unique_ptr L1L1TrackJetProducer(new l1t::TkJetWordCollection); - - // For TTStubs - const TrackerTopology &tTopo = iSetup.getData(tTopoToken_); - - edm::Handle>> TTTrackHandle; - iEvent.getByToken(trackToken_, TTTrackHandle); - vector>::const_iterator iterL1Track; - - edm::Handle PVtx; - iEvent.getByToken(PVtxToken_, PVtx); - float PVz = (PVtx->at(0)).z0(); - - L1TrkPtrs_.clear(); - zBinCount_.clear(); - tdtrk_.clear(); - unsigned int this_l1track = 0; - for (iterL1Track = TTTrackHandle->begin(); iterL1Track != TTTrackHandle->end(); iterL1Track++) { - edm::Ptr trkPtr(TTTrackHandle, this_l1track); - this_l1track++; - float trk_pt = trkPtr->momentum().perp(); - int trk_nstubs = (int)trkPtr->getStubRefs().size(); - float trk_chi2dof = trkPtr->chi2Red(); - float trk_d0 = trkPtr->d0(); - float trk_bendchi2 = trkPtr->stubPtConsistency(); - float trk_z0 = trkPtr->z0(); - - int trk_nPS = 0; - for (int istub = 0; istub < trk_nstubs; istub++) { // loop over the stubs - DetId detId(trkPtr->getStubRefs().at(istub)->getDetId()); - if (detId.det() == DetId::Detector::Tracker) { - if ((detId.subdetId() == StripSubdetector::TOB && tTopo.tobLayer(detId) <= 3) || - (detId.subdetId() == StripSubdetector::TID && tTopo.tidRing(detId) <= 9)) - trk_nPS++; - } - } - - if (trk_nPS < trkNPSStubMin_) - continue; - if (!trackQualityCuts(trk_nstubs, trk_chi2dof, trk_bendchi2)) - continue; - if (std::abs(trk_z0 - PVz) > MaxDzTrackPV) - continue; - if (std::abs(trk_z0) > trkZMax_) - continue; - if (std::abs(trkPtr->momentum().eta()) > trkEtaMax_) - continue; - if (trk_pt < trkPtMin_) - continue; - L1TrkPtrs_.push_back(trkPtr); - zBinCount_.push_back(0); - - if ((std::abs(trk_d0) > d0CutNStubs5_ && trk_nstubs >= 5 && d0CutNStubs5_ >= 0) || - (trk_nstubs == 4 && std::abs(trk_d0) > d0CutNStubs4_ && d0CutNStubs4_ >= 0)) - tdtrk_.push_back(1); //displaced track - else - tdtrk_.push_back(0); // not displaced track - } - - if (!L1TrkPtrs_.empty()) { - TrackJetEmulationMaxZBin mzb; - - L2_cluster(L1TrkPtrs_, mzb); - if (mzb.clusters != nullptr) { - for (int j = 0; j < mzb.nclust; ++j) { - //FILL Two Layer Jets for Jet Collection - if (mzb.clusters[j].pTtot < convert::makePtFromFloat(trkPtMin_)) - continue; //protects against reading bad memory - if (mzb.clusters[j].ntracks < 1) - continue; - if (mzb.clusters[j].ntracks > 5000) - continue; - l1t::TkJetWord::glbeta_t jetEta = mzb.clusters[j].eta * convert::ETA_LSB_POW; - l1t::TkJetWord::glbphi_t jetPhi = mzb.clusters[j].phi * convert::PHI_LSB_POW; - l1t::TkJetWord::z0_t jetZ0 = mzb.zbincenter * convert::Z0_LSB_POW; - l1t::TkJetWord::pt_t jetPt = mzb.clusters[j].pTtot; - l1t::TkJetWord::nt_t totalntracks_ = mzb.clusters[j].ntracks; - l1t::TkJetWord::nx_t totalxtracks_ = mzb.clusters[j].nxtracks; - l1t::TkJetWord::tkjetunassigned_t unassigned_ = 0; - - l1t::TkJetWord trkJet(jetPt, jetEta, jetPhi, jetZ0, totalntracks_, totalxtracks_, unassigned_); - //trkJet.setDispCounters(DispCounters); - L1L1TrackJetProducer->push_back(trkJet); - } - } else if (mzb.clusters == nullptr) { - edm::LogWarning("L1TrackJetEmulationProducer") << "mzb.clusters Not Assigned!\n"; - } - - if (displaced_) - iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJetsExtended"); - else - iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJets"); - delete[] mzb.clusters; - } else if (L1TrkPtrs_.empty()) { - edm::LogWarning("L1TrackJetEmulationProducer") << "L1TrkPtrs Not Assigned!\n"; - if (displaced_) - iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJetsExtended"); - else - iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJets"); - } -} - -void L1TrackJetEmulationProducer::L2_cluster(vector> L1TrkPtrs_, TrackJetEmulationMaxZBin &mzb) { - enum TrackBitWidths { - kEtaSize = 16, // Width of z-position (40cm / 0.1) - kEtaMagSize = 3, // Width of z-position magnitude (signed) - kPtSize = 14, // Width of pt - kPtMagSize = 9, // Width of pt magnitude (unsigned) - kPhiSize = 12, - kPhiMagSize = 1, - }; - - const int nz = zBins_ + 1; - TrackJetEmulationMaxZBin all_zBins[nz]; - static TrackJetEmulationMaxZBin mzbtemp; - for (int z = 0; z < nz; ++z) - all_zBins[z] = mzbtemp; - - z0_intern zmin = convert::makeZ0(-1.0 * trkZMax_); - z0_intern zmax = zmin + 2 * zStep_; - - TrackJetEmulationEtaPhiBin epbins[phiBins_][etaBins_]; // create grid of phiBins - glbphi_intern phi = convert::makeGlbPhi(-1.0 * M_PI); - glbeta_intern eta; - glbeta_intern etamin, etamax, phimin, phimax; - for (int i = 0; i < phiBins_; ++i) { - eta = convert::makeGlbEta(-1.0 * trkEtaMax_); - for (int j = 0; j < etaBins_; ++j) { - phimin = phi; - phimax = phi + phiStep_; - etamin = eta; - eta = eta + etaStep_; - etamax = eta; - epbins[i][j].phi = (phimin + phimax) / 2; - epbins[i][j].eta = (etamin + etamax) / 2; - epbins[i][j].pTtot = 0; - epbins[i][j].ntracks = 0; - epbins[i][j].nxtracks = 0; - } // for each etabin - phi = phi + phiStep_; - } // for each phibin (finished creating epbins) - - mzb = all_zBins[0]; - mzb.ht = 0; - int ntracks = L1TrkPtrs_.size(); - // uninitalized arrays - TrackJetEmulationEtaPhiBin *L1clusters[phiBins_]; - TrackJetEmulationEtaPhiBin L2cluster[ntracks]; - - for (int zbin = 0; zbin < zBins_; ++zbin) { - for (int i = 0; i < phiBins_; ++i) { //First initialize pT, numtracks, used to 0 (or false) - for (int j = 0; j < etaBins_; ++j) { - epbins[i][j].pTtot = 0; - epbins[i][j].used = false; - epbins[i][j].ntracks = 0; - epbins[i][j].nxtracks = 0; - } //for each etabin - L1clusters[i] = epbins[i]; - } //for each phibin - - for (unsigned int k = 0; k < L1TrkPtrs_.size(); ++k) { - ap_ufixed inputTrkPt = 0; - inputTrkPt.V = L1TrkPtrs_[k]->getTrackWord()(TTTrack_TrackWord::TrackBitLocations::kRinvMSB - 1, - TTTrack_TrackWord::TrackBitLocations::kRinvLSB); - pt_intern trkpt = inputTrkPt; - - ap_fixed trketainput = 0; - trketainput.V = L1TrkPtrs_[k]->getTrackWord()(TTTrack_TrackWord::TrackBitLocations::kTanlMSB, - TTTrack_TrackWord::TrackBitLocations::kTanlLSB); - ap_ufixed<64 + ETA_EXTRABITS, 32 + ETA_EXTRABITS> eta_conv = - 1.0 / convert::ETA_LSB; //conversion factor from input eta format to internal format - glbeta_intern trketa = eta_conv * trketainput; - - int Sector = L1TrkPtrs_[k]->phiSector(); - glbphi_intern trkphiSector = 0; - if (Sector < 5) { - trkphiSector = Sector * convert::makeGlbPhi(2.0 * M_PI / 9.0); - } else { - trkphiSector = - convert::makeGlbPhi(-1.0 * M_PI + M_PI / 9.0) + (Sector - 5) * convert::makeGlbPhi(2.0 * M_PI / 9.0); - } - - ap_fixed trkphiinput = 0; - trkphiinput.V = L1TrkPtrs_[k]->getTrackWord()(TTTrack_TrackWord::TrackBitLocations::kPhiMSB, - TTTrack_TrackWord::TrackBitLocations::kPhiLSB); - ap_ufixed<64 + PHI_EXTRABITS, 32 + PHI_EXTRABITS> phi_conv = - TTTrack_TrackWord::stepPhi0 / convert::PHI_LSB * - (1 << (TrackBitWidths::kPhiSize - TrackBitWidths::kPhiMagSize)); - //phi_conv is a conversion factor from input phi format to the internal format - glbeta_intern localphi = phi_conv * trkphiinput; - glbeta_intern trkphi = localphi + trkphiSector; - - ap_int inputTrkZ0 = L1TrkPtrs_[k]->getTrackWord()( - TTTrack_TrackWord::TrackBitLocations::kZ0MSB, TTTrack_TrackWord::TrackBitLocations::kZ0LSB); - ap_ufixed<32 + Z0_EXTRABITS, 1 + Z0_EXTRABITS> z0_conv = - TTTrack_TrackWord::stepZ0 / convert::Z0_LSB; //conversion factor from input z format to internal format - z0_intern trkZ = z0_conv * inputTrkZ0; - - for (int i = 0; i < phiBins_; ++i) { - for (int j = 0; j < etaBins_; ++j) { - L2cluster[k] = epbins[i][j]; - if ((zmin <= trkZ && zmax >= trkZ) && - ((epbins[i][j].eta - etaStep_ / 2 < trketa && epbins[i][j].eta + etaStep_ / 2 >= trketa) && - epbins[i][j].phi - phiStep_ / 2 < trkphi && epbins[i][j].phi + phiStep_ / 2 >= trkphi && - (zBinCount_[k] != 2))) { - zBinCount_.at(k) = zBinCount_.at(k) + 1; - - if (trkpt < convert::makePtFromFloat(trkPtMax_)) - epbins[i][j].pTtot += trkpt; - else - epbins[i][j].pTtot += convert::makePtFromFloat(trkPtMax_); - ++epbins[i][j].ntracks; - //x-bit is currently not used in firmware, so we leave nxtracks = 0 for now - } // if right bin - } // for each phibin: j loop - } // for each phibin: i loop - } // end loop over tracks - - for (int phislice = 0; phislice < phiBins_; ++phislice) { - L1clusters[phislice] = L1_cluster(epbins[phislice]); - for (int ind = 0; L1clusters[phislice][ind].pTtot != 0; ++ind) { - L1clusters[phislice][ind].used = false; - } - } - - //Create clusters array to hold output cluster data for Layer2; can't have more clusters than tracks. - //Find eta-phibin with maxpT, make center of cluster, add neighbors if not already used. - pt_intern hipT = 0; - int nclust = 0; - int phibin = 0; - int imax = -1; - int index1; //index of clusters array for each phislice - pt_intern E1 = 0; - pt_intern E0 = 0; - pt_intern E2 = 0; - l1t::TkJetWord::nt_t ntrk1, ntrk2; - l1t::TkJetWord::nx_t nxtrk1, nxtrk2; - int used1, used2, used3, used4; - - for (phibin = 0; phibin < phiBins_; ++phibin) { //Find eta-phibin with highest pT - while (true) { - hipT = 0; - for (index1 = 0; L1clusters[phibin][index1].pTtot > 0; ++index1) { - if (!L1clusters[phibin][index1].used && L1clusters[phibin][index1].pTtot >= hipT) { - hipT = L1clusters[phibin][index1].pTtot; - imax = index1; - } - } // for each index within the phibin - - if (hipT == 0) - break; // If highest pT is 0, all bins are used - E0 = hipT; // E0 is pT of first phibin of the cluster - E1 = 0; - E2 = 0; - ntrk1 = 0; - ntrk2 = 0; - nxtrk1 = 0; - nxtrk2 = 0; - L2cluster[nclust] = L1clusters[phibin][imax]; - - L1clusters[phibin][imax].used = true; - // Add pT of upper neighbor - // E1 is pT of the middle phibin (should be highest pT) - if (phibin != phiBins_ - 1) { - used1 = -1; - used2 = -1; - for (index1 = 0; L1clusters[phibin + 1][index1].pTtot != 0; ++index1) { - if (L1clusters[phibin + 1][index1].used) - continue; - if (L1clusters[phibin + 1][index1].eta - L1clusters[phibin][imax].eta <= 3 * etaStep_ / 2 && - L1clusters[phibin][imax].eta - L1clusters[phibin + 1][index1].eta <= 3 * etaStep_ / 2) { - E1 += L1clusters[phibin + 1][index1].pTtot; - ntrk1 += L1clusters[phibin + 1][index1].ntracks; - nxtrk1 += L1clusters[phibin + 1][index1].nxtracks; - if (used1 < 0) - used1 = index1; - else - used2 = index1; - } // if cluster is within one phibin - } // for each cluster in above phibin - - if (E1 < E0) { // if E1 isn't higher, E0 and E1 are their own cluster - L2cluster[nclust].pTtot += E1; - L2cluster[nclust].ntracks += ntrk1; - L2cluster[nclust].nxtracks += nxtrk1; - if (used1 >= 0) - L1clusters[phibin + 1][used1].used = true; - if (used2 >= 0) - L1clusters[phibin + 1][used2].used = true; - nclust++; - continue; - } - - if (phibin != phiBins_ - 2) { // E2 will be the pT of the third phibin (should be lower than E1) - used3 = -1; - used4 = -1; - for (index1 = 0; L1clusters[phibin + 2][index1].pTtot != 0; ++index1) { - if (L1clusters[phibin + 2][index1].used) - continue; - if (L1clusters[phibin + 2][index1].eta - L1clusters[phibin][imax].eta <= 3 * etaStep_ / 2 && - L1clusters[phibin][imax].eta - L1clusters[phibin + 2][index1].eta <= 3 * etaStep_ / 2) { - E2 += L1clusters[phibin + 2][index1].pTtot; - ntrk2 += L1clusters[phibin + 2][index1].ntracks; - nxtrk2 += L1clusters[phibin + 2][index1].nxtracks; - if (used3 < 0) - used3 = index1; - else - used4 = index1; - } - } - // if indeed E2 < E1, add E1 and E2 to E0, they're all a cluster together - // otherwise, E0 is its own cluster - if (E2 < E1) { - L2cluster[nclust].pTtot += E1 + E2; - L2cluster[nclust].ntracks += ntrk1 + ntrk2; - L2cluster[nclust].nxtracks += nxtrk1 + nxtrk2; - L2cluster[nclust].phi = L1clusters[phibin + 1][used1].phi; - if (used1 >= 0) - L1clusters[phibin + 1][used1].used = true; - if (used2 >= 0) - L1clusters[phibin + 1][used2].used = true; - if (used3 >= 0) - L1clusters[phibin + 2][used3].used = true; - if (used4 >= 0) - L1clusters[phibin + 2][used4].used = true; - } - nclust++; - continue; - } // end Not phiBins-2 - else { - L2cluster[nclust].pTtot += E1; - L2cluster[nclust].ntracks += ntrk1; - L2cluster[nclust].nxtracks += nxtrk1; - L2cluster[nclust].phi = L1clusters[phibin + 1][used1].phi; - if (used1 >= 0) - L1clusters[phibin + 1][used1].used = true; - if (used2 >= 0) - L1clusters[phibin + 1][used2].used = true; - nclust++; - continue; - } - } //End not last phibin(23) - else { //if it is phibin 23 - L1clusters[phibin][imax].used = true; - nclust++; - } - } // while hipT not 0 - } // for each phibin - - for (phibin = 0; phibin < phiBins_; ++phibin) - delete[] L1clusters[phibin]; - - // Now merge clusters, if necessary - for (int m = 0; m < nclust - 1; ++m) { - for (int n = m + 1; n < nclust; ++n) { - if (L2cluster[n].eta == L2cluster[m].eta && - ((L2cluster[n].phi - L2cluster[m].phi < 3 * phiStep_ / 2 && - L2cluster[m].phi - L2cluster[n].phi < 3 * phiStep_ / 2) || - (L2cluster[n].phi - L2cluster[m].phi > 2 * convert::makeGlbPhi(M_PI) - phiStep_ || - L2cluster[m].phi - L2cluster[n].phi > 2 * convert::makeGlbPhi(M_PI) - phiStep_))) { - if (L2cluster[n].pTtot > L2cluster[m].pTtot) { - L2cluster[m].phi = L2cluster[n].phi; - } - L2cluster[m].pTtot += L2cluster[n].pTtot; - L2cluster[m].ntracks += L2cluster[n].ntracks; - L2cluster[m].nxtracks += L2cluster[n].nxtracks; - for (int m1 = n; m1 < nclust - 1; ++m1) { - L2cluster[m1] = L2cluster[m1 + 1]; - } - nclust--; - m = -1; - break; //????? - } // end if clusters neighbor in eta - } - } // end for (m) loop - // sum up all pTs in this zbin to find ht - - pt_intern ht = 0; - for (int k = 0; k < nclust; ++k) { - if (L2cluster[k].pTtot > convert::makePtFromFloat(lowpTJetMinpT_) && - L2cluster[k].ntracks < lowpTJetMinTrackMultiplicity_) - continue; - if (L2cluster[k].pTtot > convert::makePtFromFloat(highpTJetMinpT_) && - L2cluster[k].ntracks < highpTJetMinTrackMultiplicity_) - continue; - if (L2cluster[k].pTtot > convert::makePtFromFloat(minTrkJetpT_)) - ht += L2cluster[k].pTtot; - } - // if ht is larger than previous max, this is the new vertex zbin - all_zBins[zbin].znum = zbin; - all_zBins[zbin].clusters = new TrackJetEmulationEtaPhiBin[nclust]; - all_zBins[zbin].nclust = nclust; - all_zBins[zbin].zbincenter = (zmin + zmax) / 2; - for (int k = 0; k < nclust; ++k) { - all_zBins[zbin].clusters[k].phi = L2cluster[k].phi; - all_zBins[zbin].clusters[k].eta = L2cluster[k].eta; - all_zBins[zbin].clusters[k].pTtot = L2cluster[k].pTtot; - all_zBins[zbin].clusters[k].ntracks = L2cluster[k].ntracks; - all_zBins[zbin].clusters[k].nxtracks = L2cluster[k].nxtracks; - } - all_zBins[zbin].ht = ht; - if (ht >= mzb.ht) { - mzb = all_zBins[zbin]; - mzb.zbincenter = (zmin + zmax) / 2; - } - // Prepare for next zbin! - zmin = zmin + zStep_; - zmax = zmax + zStep_; - } // for each zbin - - for (int zbin = 0; zbin < zBins_; ++zbin) { - if (zbin == mzb.znum) { - continue; - } - delete[] all_zBins[zbin].clusters; - } -} - -TrackJetEmulationEtaPhiBin *L1TrackJetEmulationProducer::L1_cluster(TrackJetEmulationEtaPhiBin *phislice) { - TrackJetEmulationEtaPhiBin *clusters = new TrackJetEmulationEtaPhiBin[etaBins_ / 2]; - for (int etabin = 0; etabin < etaBins_ / 2; ++etabin) { - clusters[etabin].pTtot = 0; - clusters[etabin].ntracks = 0; - clusters[etabin].nxtracks = 0; - clusters[etabin].phi = 0; - clusters[etabin].eta = 0; - clusters[etabin].used = false; - } - - if (clusters == nullptr) - edm::LogWarning("L1TrackJetEmulationProducer") << "Clusters memory not assigned!\n"; - - // Find eta-phibin with maxpT, make center of cluster, add neighbors if not already used - pt_intern my_pt, left_pt, right_pt, right2pt; - int nclust = 0; - right2pt = 0; - for (int etabin = 0; etabin < etaBins_; ++etabin) { - // assign values for my pT and neighbors' pT - if (phislice[etabin].used) - continue; - my_pt = phislice[etabin].pTtot; - if (etabin > 0 && !phislice[etabin - 1].used) { - left_pt = phislice[etabin - 1].pTtot; - } else - left_pt = 0; - if (etabin < etaBins_ - 1 && !phislice[etabin + 1].used) { - right_pt = phislice[etabin + 1].pTtot; - if (etabin < etaBins_ - 2 && !phislice[etabin + 2].used) { - right2pt = phislice[etabin + 2].pTtot; - } else - right2pt = 0; - } else - right_pt = 0; - - // if I'm not a cluster, move on - if (my_pt < left_pt || my_pt <= right_pt) { - // if unused pT in the left neighbor, spit it out as a cluster - if (left_pt > 0) { - clusters[nclust] = phislice[etabin - 1]; - phislice[etabin - 1].used = true; - nclust++; - } - continue; - } - - // I guess I'm a cluster-- should I use my right neighbor? - // Note: left neighbor will definitely be used because if it - // didn't belong to me it would have been used already - clusters[nclust] = phislice[etabin]; - phislice[etabin].used = true; - if (left_pt > 0) { - clusters[nclust].pTtot += left_pt; - clusters[nclust].ntracks += phislice[etabin - 1].ntracks; - clusters[nclust].nxtracks += phislice[etabin - 1].nxtracks; - } - if (my_pt >= right2pt && right_pt > 0) { - clusters[nclust].pTtot += right_pt; - clusters[nclust].ntracks += phislice[etabin + 1].ntracks; - clusters[nclust].nxtracks += phislice[etabin + 1].nxtracks; - phislice[etabin + 1].used = true; - } - nclust++; - } // for each etabin - - // Now merge clusters, if necessary - for (int m = 0; m < nclust - 1; ++m) { - if (clusters[m + 1].eta - clusters[m].eta < 3 * etaStep_ / 2 && - clusters[m].eta - clusters[m + 1].eta < 3 * etaStep_ / 2) { - if (clusters[m + 1].pTtot > clusters[m].pTtot) { - clusters[m].eta = clusters[m + 1].eta; - } - clusters[m].pTtot += clusters[m + 1].pTtot; - clusters[m].ntracks += clusters[m + 1].ntracks; // Previous version didn't add tracks when merging - clusters[m].nxtracks += clusters[m + 1].nxtracks; - for (int m1 = m + 1; m1 < nclust - 1; ++m1) - clusters[m1] = clusters[m1 + 1]; - nclust--; - m = -1; - } // end if clusters neighbor in eta - } // end for (m) loop - - for (int i = nclust; i < etaBins_ / 2; ++i) // zero out remaining unused clusters - clusters[i].pTtot = 0; - return clusters; -} - -void L1TrackJetEmulationProducer::beginStream(StreamID) {} - -void L1TrackJetEmulationProducer::endStream() {} - -bool L1TrackJetEmulationProducer::trackQualityCuts(int trk_nstub, float trk_chi2, float trk_bendchi2) { - bool PassQuality = false; - if (trk_bendchi2 < trkBendChi2Max_ && trk_chi2 < trkChi2dofMax_ && trk_nstub >= 4 && !displaced_) - PassQuality = true; - if (displaced_ && trk_bendchi2 < nStubs4DisplacedBend_ && trk_chi2 < nStubs4DisplacedChi2_ && trk_nstub == 4) - PassQuality = true; - if (displaced_ && trk_bendchi2 < nStubs5DisplacedBend_ && trk_chi2 < nStubs5DisplacedChi2_ && trk_nstub > 4) - PassQuality = true; - return PassQuality; -} - -void L1TrackJetEmulationProducer::fillDescriptions(ConfigurationDescriptions &descriptions) { - //The following says we do not know what parameters are allowed so do no validation - // Please change this to state exactly what you do use, even if it is no parameters - ParameterSetDescription desc; - desc.setUnknown(); - descriptions.addDefault(desc); -} - -//define this as a plug-in -DEFINE_FWK_MODULE(L1TrackJetEmulationProducer); diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetEmulatorProducer.cc b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetEmulatorProducer.cc new file mode 100644 index 0000000000000..96ce847d7efa5 --- /dev/null +++ b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetEmulatorProducer.cc @@ -0,0 +1,456 @@ +// Original Author: Samuel Edwin Leigh, Tyler Wu, +// Rutgers, the State University of New Jersey +// +// Rewritting/Improvements: George Karathanasis, +// georgios.karathanasis@cern.ch, CU Boulder +// +// Created: Wed, 01 Aug 2018 14:01:41 GMT +// Latest update: Nov 2022 (by GK) +// +// Track jets are clustered in a two-layer process, first by clustering in phi, +// then by clustering in eta. The code proceeds as following: putting all tracks// in a grid of eta vs phi space, and then cluster them. Finally we merge the cl +// usters when needed. The code is improved to use the same module between emula +// tion and simulation was also improved, with bug fixes and being faster. +// Introduction to object (p10-13): +// https://indico.cern.ch/event/791517/contributions/3341650/attachments/1818736/2973771/TrackBasedAlgos_L1TMadrid_MacDonald.pdf +// New and improved version: https://indico.cern.ch/event/1203796/contributions/5073056/attachments/2519806/4333006/trackjet_emu.pdf + +// L1T include files +#include "DataFormats/L1TCorrelator/interface/TkJet.h" +#include "DataFormats/L1TCorrelator/interface/TkJetFwd.h" +#include "DataFormats/L1TrackTrigger/interface/TTTypes.h" +#include "DataFormats/L1TrackTrigger/interface/TTTrack.h" +#include "DataFormats/L1TrackTrigger/interface/TTTrack_TrackWord.h" +#include "DataFormats/L1Trigger/interface/TkJetWord.h" +#include "DataFormats/L1Trigger/interface/VertexWord.h" + +// system include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "DataFormats/Common/interface/Ref.h" + +//own headers +#include "L1TrackJetClustering.h" + +//general +#include + +using namespace std; +using namespace edm; +using namespace l1t; +using namespace l1ttrackjet; + +class L1TrackJetEmulatorProducer : public stream::EDProducer<> { +public: + explicit L1TrackJetEmulatorProducer(const ParameterSet &); + ~L1TrackJetEmulatorProducer() override = default; + typedef TTTrack L1TTTrackType; + typedef vector L1TTTrackCollectionType; + static void fillDescriptions(ConfigurationDescriptions &descriptions); + +private: + void produce(Event &, const EventSetup &) override; + + // ----------member data --------------------------- + + std::vector> L1TrkPtrs_; + vector tdtrk_; + const float trkZMax_; + const float trkPtMax_; + const float trkPtMin_; + const float trkEtaMax_; + const float nStubs4PromptChi2_; + const float nStubs5PromptChi2_; + const float nStubs4PromptBend_; + const float nStubs5PromptBend_; + const int trkNPSStubMin_; + const int lowpTJetMinTrackMultiplicity_; + const float lowpTJetThreshold_; + const int highpTJetMinTrackMultiplicity_; + const float highpTJetThreshold_; + const int zBins_; + const int etaBins_; + const int phiBins_; + const double minTrkJetpT_; + const bool displaced_; + const float d0CutNStubs4_; + const float d0CutNStubs5_; + const float nStubs4DisplacedChi2_; + const float nStubs5DisplacedChi2_; + const float nStubs4DisplacedBend_; + const float nStubs5DisplacedBend_; + const int nDisplacedTracks_; + const float dzPVTrk_; + + float PVz; + float zStep_; + glbeta_intern etaStep_; + glbphi_intern phiStep_; + + TTTrack_TrackWord trackword; + + edm::ESGetToken tTopoToken_; + const EDGetTokenT>> trackToken_; + const EDGetTokenT PVtxToken_; +}; + +//constructor +L1TrackJetEmulatorProducer::L1TrackJetEmulatorProducer(const ParameterSet &iConfig) + : trkZMax_(iConfig.getParameter("trk_zMax")), + trkPtMax_(iConfig.getParameter("trk_ptMax")), + trkPtMin_(iConfig.getParameter("trk_ptMin")), + trkEtaMax_(iConfig.getParameter("trk_etaMax")), + nStubs4PromptChi2_(iConfig.getParameter("nStubs4PromptChi2")), + nStubs5PromptChi2_(iConfig.getParameter("nStubs5PromptChi2")), + nStubs4PromptBend_(iConfig.getParameter("nStubs4PromptBend")), + nStubs5PromptBend_(iConfig.getParameter("nStubs5PromptBend")), + trkNPSStubMin_(iConfig.getParameter("trk_nPSStubMin")), + lowpTJetMinTrackMultiplicity_(iConfig.getParameter("lowpTJetMinTrackMultiplicity")), + lowpTJetThreshold_(iConfig.getParameter("lowpTJetThreshold")), + highpTJetMinTrackMultiplicity_(iConfig.getParameter("highpTJetMinTrackMultiplicity")), + highpTJetThreshold_(iConfig.getParameter("highpTJetThreshold")), + zBins_(iConfig.getParameter("zBins")), + etaBins_(iConfig.getParameter("etaBins")), + phiBins_(iConfig.getParameter("phiBins")), + minTrkJetpT_(iConfig.getParameter("minTrkJetpT")), + displaced_(iConfig.getParameter("displaced")), + d0CutNStubs4_(iConfig.getParameter("d0_cutNStubs4")), + d0CutNStubs5_(iConfig.getParameter("d0_cutNStubs5")), + nStubs4DisplacedChi2_(iConfig.getParameter("nStubs4DisplacedChi2")), + nStubs5DisplacedChi2_(iConfig.getParameter("nStubs5DisplacedChi2")), + nStubs4DisplacedBend_(iConfig.getParameter("nStubs4DisplacedBend")), + nStubs5DisplacedBend_(iConfig.getParameter("nStubs5DisplacedBend")), + nDisplacedTracks_(iConfig.getParameter("nDisplacedTracks")), + dzPVTrk_(iConfig.getParameter("MaxDzTrackPV")), + tTopoToken_(esConsumes(edm::ESInputTag("", ""))), + trackToken_(consumes>>(iConfig.getParameter("L1TrackInputTag"))), + PVtxToken_(consumes(iConfig.getParameter("L1PVertexInputTag"))) { + zStep_ = 2.0 * trkZMax_ / (zBins_ + 1); // added +1 in denom + etaStep_ = glbeta_intern(2.0 * trkEtaMax_ / etaBins_); //etaStep is the width of an etabin + phiStep_ = DoubleToBit(2.0 * (M_PI) / phiBins_, + TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, + TTTrack_TrackWord::stepPhi0); ///phiStep is the width of a phibin + + if (displaced_) + produces("L1TrackJetsExtended"); + else + produces("L1TrackJets"); +} + +void L1TrackJetEmulatorProducer::produce(Event &iEvent, const EventSetup &iSetup) { + unique_ptr L1TrackJetContainer(new l1t::TkJetWordCollection); + // Read inputs + const TrackerTopology &tTopo = iSetup.getData(tTopoToken_); + + edm::Handle>> TTTrackHandle; + iEvent.getByToken(trackToken_, TTTrackHandle); + + edm::Handle PVtx; + iEvent.getByToken(PVtxToken_, PVtx); + float PVz = (PVtx->at(0)).z0(); + + L1TrkPtrs_.clear(); + tdtrk_.clear(); + // track selection + for (unsigned int this_l1track = 0; this_l1track < TTTrackHandle->size(); this_l1track++) { + edm::Ptr trkPtr(TTTrackHandle, this_l1track); + float trk_pt = trkPtr->momentum().perp(); + int trk_nstubs = (int)trkPtr->getStubRefs().size(); + float trk_chi2dof = trkPtr->chi2Red(); + float trk_bendchi2 = trkPtr->stubPtConsistency(); + int trk_nPS = 0; + for (int istub = 0; istub < trk_nstubs; istub++) { + DetId detId(trkPtr->getStubRefs().at(istub)->getDetId()); + if (detId.det() == DetId::Detector::Tracker) { + if ((detId.subdetId() == StripSubdetector::TOB && tTopo.tobLayer(detId) <= 3) || + (detId.subdetId() == StripSubdetector::TID && tTopo.tidRing(detId) <= 9)) + trk_nPS++; + } + } + // selection tracks - supposed to happen on seperate module (kept for legacy/debug reasons) + if (trk_nPS < trkNPSStubMin_) + continue; + if (!TrackQualitySelection(trk_nstubs, + trk_chi2dof, + trk_bendchi2, + nStubs4PromptBend_, + nStubs5PromptBend_, + nStubs4PromptChi2_, + nStubs5PromptChi2_, + nStubs4DisplacedBend_, + nStubs5DisplacedBend_, + nStubs4DisplacedChi2_, + nStubs5DisplacedChi2_, + displaced_)) + continue; + if (std::abs(PVz - trkPtr->z0()) > dzPVTrk_ && dzPVTrk_ > 0) + continue; + if (std::abs(trkPtr->z0()) > trkZMax_) + continue; + if (std::abs(trkPtr->momentum().eta()) > trkEtaMax_) + continue; + if (trk_pt < trkPtMin_) + continue; + L1TrkPtrs_.push_back(trkPtr); + } + + // if no tracks pass selection return empty containers + if (L1TrkPtrs_.empty()) { + if (displaced_) + iEvent.put(std::move(L1TrackJetContainer), "L1TrackJetsExtended"); + else + iEvent.put(std::move(L1TrackJetContainer), "L1TrackJets"); + return; + } + + TrackJetEmulationMaxZBin mzb; + mzb.znum = 0; + mzb.nclust = 0; + mzb.ht = 0; + + TrackJetEmulationEtaPhiBin epbins_default[phiBins_][etaBins_]; // create grid of phiBins + glbphi_intern phi = DoubleToBit( + -1.0 * M_PI, TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, TTTrack_TrackWord::stepPhi0); + for (int i = 0; i < phiBins_; ++i) { + glbeta_intern eta = -1 * trkEtaMax_; + for (int j = 0; j < etaBins_; ++j) { + epbins_default[i][j].phi = (phi + (phi + phiStep_)) / 2; // phi coord of bin center + epbins_default[i][j].eta = (eta + (eta + etaStep_)) / 2; // eta coord of bin center + eta += etaStep_; + } // for each etabin + phi += phiStep_; + } // for each phibin (finished creating epbins) + + // bins for z if multibin option is selected + std::vector zmins, zmaxs; + for (int zbin = 0; zbin < zBins_; zbin++) { + zmins.push_back(DoubleToBit( + -1.0 * trkZMax_ + zStep_ * zbin, TTTrack_TrackWord::TrackBitWidths::kZ0Size, TTTrack_TrackWord::stepZ0)); + zmaxs.push_back(DoubleToBit(-1.0 * trkZMax_ + zStep_ * zbin + 2 * zStep_, + TTTrack_TrackWord::TrackBitWidths::kZ0Size, + TTTrack_TrackWord::stepZ0)); + } + + // create vectors that hold clusters + std::vector> L1clusters; + L1clusters.reserve(phiBins_); + std::vector L2clusters; + + // default (empty) grid + for (int i = 0; i < phiBins_; ++i) { + for (int j = 0; j < etaBins_; ++j) { + epbins_default[i][j].pTtot = 0; + epbins_default[i][j].used = false; + epbins_default[i][j].ntracks = 0; + epbins_default[i][j].nxtracks = 0; + epbins_default[i][j].trackidx.clear(); + } + } + + // logic: loop over z bins find tracks in this bin and arrange them in a 2D eta-phi matrix + for (unsigned int zbin = 0; zbin < zmins.size(); ++zbin) { + // initialize matrices for every z bin + z0_intern zmin = zmins[zbin]; + z0_intern zmax = zmaxs[zbin]; + TrackJetEmulationEtaPhiBin epbins[phiBins_][etaBins_]; + std::copy(&epbins_default[0][0], &epbins_default[0][0] + phiBins_ * etaBins_, &epbins[0][0]); + + //clear containers + L1clusters.clear(); + L2clusters.clear(); + + // fill grid + for (unsigned int k = 0; k < L1TrkPtrs_.size(); ++k) { + //// conversions + //-z0 + z0_intern trkZ = L1TrkPtrs_[k]->getZ0Word(); + + if (zmax < trkZ) + continue; + if (zmin > trkZ) + continue; + if (zbin == 0 && zmin == trkZ) + continue; + + //-pt + ap_uint ptEmulationBits = L1TrkPtrs_[k]->getRinvWord(); + pt_intern trkpt; + trkpt.V = ptEmulationBits.range(); + + //-eta + TTTrack_TrackWord::tanl_t etaEmulationBits = L1TrkPtrs_[k]->getTanlWord(); + glbeta_intern trketa; + trketa.V = etaEmulationBits.range(); + + //-phi + int Sector = L1TrkPtrs_[k]->phiSector(); + double sector_phi_value = 0; + if (Sector < 5) { + sector_phi_value = 2.0 * M_PI * Sector / 9.0; + } else { + sector_phi_value = (-1.0 * M_PI + M_PI / 9.0 + (Sector - 5) * 2.0 * M_PI / 9.0); + } + + glbphi_intern trkphiSector = DoubleToBit(sector_phi_value, + TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, + TTTrack_TrackWord::stepPhi0); + glbphi_intern local_phi = 0; + local_phi.V = L1TrkPtrs_[k]->getPhiWord(); + glbphi_intern local_phi2 = + DoubleToBit(BitToDouble(local_phi, TTTrack_TrackWord::TrackBitWidths::kPhiSize, TTTrack_TrackWord::stepPhi0), + TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, + TTTrack_TrackWord::stepPhi0); + glbphi_intern trkphi = local_phi2 + trkphiSector; + + //-d0 + d0_intern abs_trkD0 = L1TrkPtrs_[k]->getD0Word(); + + //-nstub + int trk_nstubs = (int)L1TrkPtrs_[k]->getStubRefs().size(); + + // now fill the 2D grid with tracks + for (int i = 0; i < phiBins_; ++i) { + for (int j = 0; j < etaBins_; ++j) { + glbeta_intern eta_min = epbins[i][j].eta - etaStep_ / 2; //eta min + glbeta_intern eta_max = epbins[i][j].eta + etaStep_ / 2; //eta max + glbphi_intern phi_min = epbins[i][j].phi - phiStep_ / 2; //phi min + glbphi_intern phi_max = epbins[i][j].phi + phiStep_ / 2; //phi max + if ((trketa < eta_min) && j != 0) + continue; + if ((trketa > eta_max) && j != etaBins_ - 1) + continue; + if ((trkphi < phi_min) && i != 0) + continue; + if ((trkphi > phi_max) && i != phiBins_ - 1) + continue; + + if (trkpt < pt_intern(trkPtMax_)) + epbins[i][j].pTtot += trkpt; + else + epbins[i][j].pTtot += pt_intern(trkPtMax_); + if ((abs_trkD0 > + DoubleToBit(d0CutNStubs5_, TTTrack_TrackWord::TrackBitWidths::kD0Size, TTTrack_TrackWord::stepD0) && + trk_nstubs >= 5 && d0CutNStubs5_ >= 0) || + (abs_trkD0 > + DoubleToBit(d0CutNStubs4_, TTTrack_TrackWord::TrackBitWidths::kD0Size, TTTrack_TrackWord::stepD0) && + trk_nstubs == 4 && d0CutNStubs4_ >= 0)) + epbins[i][j].nxtracks += 1; + + epbins[i][j].trackidx.push_back(k); + ++epbins[i][j].ntracks; + } // for each etabin + } // for each phibin + } //end loop over tracks + + // first layer clustering - in eta using grid + for (int phibin = 0; phibin < phiBins_; ++phibin) { + L1clusters.push_back(L1_clustering( + epbins[phibin], etaBins_, etaStep_)); + } + + //second layer clustering in phi for all eta clusters + L2clusters = L2_clustering( + L1clusters, phiBins_, phiStep_, etaStep_); + + // sum all cluster pTs in this zbin to find max + pt_intern sum_pt = 0; + for (unsigned int k = 0; k < L2clusters.size(); ++k) { + if (L2clusters[k].pTtot > pt_intern(highpTJetThreshold_) && L2clusters[k].ntracks < lowpTJetMinTrackMultiplicity_) + continue; + if (L2clusters[k].pTtot > pt_intern(highpTJetThreshold_) && + L2clusters[k].ntracks < highpTJetMinTrackMultiplicity_) + continue; + + if (L2clusters[k].pTtot > pt_intern(minTrkJetpT_)) + sum_pt += L2clusters[k].pTtot; + } + if (sum_pt < mzb.ht) + continue; + // if ht is larger than previous max, this is the new vertex zbin + mzb.ht = sum_pt; + mzb.znum = zbin; + mzb.clusters = L2clusters; + mzb.nclust = L2clusters.size(); + mzb.zbincenter = (zmin + zmax) / 2.0; + } //zbin loop + + vector> L1TrackAssocJet; + for (unsigned int j = 0; j < mzb.clusters.size(); ++j) { + if (mzb.clusters[j].pTtot < pt_intern(trkPtMin_)) + continue; + + l1t::TkJetWord::glbeta_t jetEta = DoubleToBit(double(mzb.clusters[j].eta), + TkJetWord::TkJetBitWidths::kGlbEtaSize, + TkJetWord::MAX_ETA / (1 << TkJetWord::TkJetBitWidths::kGlbEtaSize)); + l1t::TkJetWord::glbphi_t jetPhi = DoubleToBit( + BitToDouble(mzb.clusters[j].phi, TTTrack_TrackWord::TrackBitWidths::kPhiSize + 4, TTTrack_TrackWord::stepPhi0), + TkJetWord::TkJetBitWidths::kGlbPhiSize, + (2. * std::abs(M_PI)) / (1 << TkJetWord::TkJetBitWidths::kGlbPhiSize)); + l1t::TkJetWord::z0_t jetZ0 = 0; + l1t::TkJetWord::pt_t jetPt = mzb.clusters[j].pTtot; + l1t::TkJetWord::nt_t total_ntracks = mzb.clusters[j].ntracks; + l1t::TkJetWord::nx_t total_disptracks = mzb.clusters[j].nxtracks; + l1t::TkJetWord::dispflag_t dispflag = 0; + l1t::TkJetWord::tkjetunassigned_t unassigned = 0; + + if (total_disptracks > nDisplacedTracks_ || total_disptracks == nDisplacedTracks_) + dispflag = 1; + L1TrackAssocJet.clear(); + for (unsigned int itrk = 0; itrk < mzb.clusters[j].trackidx.size(); itrk++) + L1TrackAssocJet.push_back(L1TrkPtrs_[mzb.clusters[j].trackidx[itrk]]); + + l1t::TkJetWord trkJet(jetPt, jetEta, jetPhi, jetZ0, total_ntracks, total_disptracks, dispflag, unassigned); + + L1TrackJetContainer->push_back(trkJet); + } + + std::sort(L1TrackJetContainer->begin(), L1TrackJetContainer->end(), [](auto &a, auto &b) { return a.pt() > b.pt(); }); + if (displaced_) + iEvent.put(std::move(L1TrackJetContainer), "L1TrackJetsExtended"); + else + iEvent.put(std::move(L1TrackJetContainer), "L1TrackJets"); +} + +void L1TrackJetEmulatorProducer::fillDescriptions(ConfigurationDescriptions &descriptions) { + //The following says we do not know what parameters are allowed so do no validation + // Please change this to state exactly what you do use, even if it is no parameters + ParameterSetDescription desc; + desc.add("L1TrackInputTag", edm::InputTag("l1tTTTracksFromTrackletEmulation", "Level1TTTracks")); + desc.add("L1PVertexInputTag", edm::InputTag("l1tVertexFinderEmulator", "l1verticesEmulation")); + desc.add("MaxDzTrackPV", 1.0); + desc.add("trk_zMax", 15.0); + desc.add("trk_ptMax", 200.0); + desc.add("trk_ptMin", 3.0); + desc.add("trk_etaMax", 2.4); + desc.add("nStubs4PromptChi2", 5.0); + desc.add("nStubs4PromptBend", 1.7); + desc.add("nStubs5PromptChi2", 2.75); + desc.add("nStubs5PromptBend", 3.5); + desc.add("trk_nPSStubMin", -1); + desc.add("minTrkJetpT", -1.0); + desc.add("etaBins", 24); + desc.add("phiBins", 27); + desc.add("zBins", 1); + desc.add("d0_cutNStubs4", -1); + desc.add("d0_cutNStubs5", -1); + desc.add("lowpTJetMinTrackMultiplicity", 2); + desc.add("lowpTJetThreshold", 50.0); + desc.add("highpTJetMinTrackMultiplicity", 3); + desc.add("highpTJetThreshold", 100.0); + desc.add("displaced", false); + desc.add("nStubs4DisplacedChi2", 5.0); + desc.add("nStubs4DisplacedBend", 1.7); + desc.add("nStubs5DisplacedChi2", 2.75); + desc.add("nStubs5DisplacedBend", 3.5); + desc.add("nDisplacedTracks", 2); + descriptions.add("l1tTrackJetsEmulator", desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(L1TrackJetEmulatorProducer); diff --git a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc index 32c40ba9d37a4..9a93ddb686ea6 100644 --- a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc +++ b/L1Trigger/L1TTrackMatch/plugins/L1TrackJetProducer.cc @@ -1,20 +1,28 @@ -// Original Author: Rishi Patel -// Modifications: George Karathanasis, georgios.karathanasis@cern.ch, CU Boulder +// Original simulation Author: Rishi Patel +// +// Rewritting/improvements: George Karathanasis, +// georgios.karathanasis@cern.ch, CU Boulder +// // Created: Wed, 01 Aug 2018 14:01:41 GMT -// Latest update: Nov 2021 (by GK) +// Latest update: Nov 2022 (by GK) // // Track jets are clustered in a two-layer process, first by clustering in phi, -// then by clustering in eta +// then by clustering in eta. The code proceeds as following: putting all tracks +// in a grid of eta vs phi space, and then cluster them. Finally we merge the cl +// usters when needed. The code is improved to use the same module between emula +// tion and simulation was also improved, with bug fixes and being faster. // Introduction to object (p10-13): // https://indico.cern.ch/event/791517/contributions/3341650/attachments/1818736/2973771/TrackBasedAlgos_L1TMadrid_MacDonald.pdf +// New and improved version: https://indico.cern.ch/event/1203796/contributions/5073056/attachments/2519806/4333006/trackjet_emu.pdf // system include files - #include "DataFormats/Common/interface/Ref.h" #include "DataFormats/L1TCorrelator/interface/TkJet.h" #include "DataFormats/L1TCorrelator/interface/TkJetFwd.h" #include "DataFormats/L1TrackTrigger/interface/TTTypes.h" #include "DataFormats/L1TrackTrigger/interface/TTTrack.h" +#include "DataFormats/L1Trigger/interface/VertexWord.h" + // user include files #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/Framework/interface/stream/EDProducer.h" @@ -24,42 +32,25 @@ #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/Utilities/interface/StreamID.h" #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" -#include "Geometry/CommonTopologies/interface/PixelGeomDetUnit.h" -#include "Geometry/CommonTopologies/interface/PixelGeomDetType.h" -#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" - -#include "DataFormats/L1Trigger/interface/Vertex.h" -#include "L1Trigger/L1TTrackMatch/interface/L1Clustering.h" - -#include "TH1D.h" -#include "TH2D.h" -#include -#include -#include -#include -#include -#include -#include +//own headers +#include "L1TrackJetClustering.h" using namespace std; using namespace edm; using namespace l1t; +using namespace l1ttrackjet; class L1TrackJetProducer : public stream::EDProducer<> { public: explicit L1TrackJetProducer(const ParameterSet &); - ~L1TrackJetProducer() override; + ~L1TrackJetProducer() override = default; typedef TTTrack L1TTTrackType; typedef vector L1TTTrackCollectionType; - static void fillDescriptions(ConfigurationDescriptions &descriptions); - bool trackQualityCuts(int trk_nstub, float trk_chi2, float trk_bendchi2); private: - void beginStream(StreamID) override; void produce(Event &, const EventSetup &) override; - void endStream() override; // ----------member data --------------------------- @@ -97,39 +88,39 @@ class L1TrackJetProducer : public stream::EDProducer<> { edm::ESGetToken tTopoToken_; const EDGetTokenT>> trackToken_; - const edm::EDGetTokenT> PVtxToken_; + const EDGetTokenT PVtxToken_; }; L1TrackJetProducer::L1TrackJetProducer(const ParameterSet &iConfig) - : trkZMax_((float)iConfig.getParameter("trk_zMax")), - trkPtMax_((float)iConfig.getParameter("trk_ptMax")), - trkPtMin_((float)iConfig.getParameter("trk_ptMin")), - trkEtaMax_((float)iConfig.getParameter("trk_etaMax")), - nStubs4PromptChi2_((float)iConfig.getParameter("nStubs4PromptChi2")), - nStubs5PromptChi2_((float)iConfig.getParameter("nStubs5PromptChi2")), - nStubs4PromptBend_((float)iConfig.getParameter("nStubs4PromptBend")), - nStubs5PromptBend_((float)iConfig.getParameter("nStubs5PromptBend")), - trkNPSStubMin_((int)iConfig.getParameter("trk_nPSStubMin")), - lowpTJetMinTrackMultiplicity_((int)iConfig.getParameter("lowpTJetMinTrackMultiplicity")), - lowpTJetThreshold_((float)iConfig.getParameter("lowpTJetThreshold")), - highpTJetMinTrackMultiplicity_((int)iConfig.getParameter("highpTJetMinTrackMultiplicity")), - highpTJetThreshold_((float)iConfig.getParameter("highpTJetThreshold")), - zBins_((int)iConfig.getParameter("zBins")), - etaBins_((int)iConfig.getParameter("etaBins")), - phiBins_((int)iConfig.getParameter("phiBins")), + : trkZMax_(iConfig.getParameter("trk_zMax")), + trkPtMax_(iConfig.getParameter("trk_ptMax")), + trkPtMin_(iConfig.getParameter("trk_ptMin")), + trkEtaMax_(iConfig.getParameter("trk_etaMax")), + nStubs4PromptChi2_(iConfig.getParameter("nStubs4PromptChi2")), + nStubs5PromptChi2_(iConfig.getParameter("nStubs5PromptChi2")), + nStubs4PromptBend_(iConfig.getParameter("nStubs4PromptBend")), + nStubs5PromptBend_(iConfig.getParameter("nStubs5PromptBend")), + trkNPSStubMin_(iConfig.getParameter("trk_nPSStubMin")), + lowpTJetMinTrackMultiplicity_(iConfig.getParameter("lowpTJetMinTrackMultiplicity")), + lowpTJetThreshold_(iConfig.getParameter("lowpTJetThreshold")), + highpTJetMinTrackMultiplicity_(iConfig.getParameter("highpTJetMinTrackMultiplicity")), + highpTJetThreshold_(iConfig.getParameter("highpTJetThreshold")), + zBins_(iConfig.getParameter("zBins")), + etaBins_(iConfig.getParameter("etaBins")), + phiBins_(iConfig.getParameter("phiBins")), minTrkJetpT_(iConfig.getParameter("minTrkJetpT")), displaced_(iConfig.getParameter("displaced")), - d0CutNStubs4_((float)iConfig.getParameter("d0_cutNStubs4")), - d0CutNStubs5_((float)iConfig.getParameter("d0_cutNStubs5")), - nStubs4DisplacedChi2_((float)iConfig.getParameter("nStubs4DisplacedChi2")), - nStubs5DisplacedChi2_((float)iConfig.getParameter("nStubs5DisplacedChi2")), - nStubs4DisplacedBend_((float)iConfig.getParameter("nStubs4DisplacedBend")), - nStubs5DisplacedBend_((float)iConfig.getParameter("nStubs5DisplacedBend")), - nDisplacedTracks_((int)iConfig.getParameter("nDisplacedTracks")), - dzPVTrk_((float)iConfig.getParameter("MaxDzTrackPV")), + d0CutNStubs4_(iConfig.getParameter("d0_cutNStubs4")), + d0CutNStubs5_(iConfig.getParameter("d0_cutNStubs5")), + nStubs4DisplacedChi2_(iConfig.getParameter("nStubs4DisplacedChi2")), + nStubs5DisplacedChi2_(iConfig.getParameter("nStubs5DisplacedChi2")), + nStubs4DisplacedBend_(iConfig.getParameter("nStubs4DisplacedBend")), + nStubs5DisplacedBend_(iConfig.getParameter("nStubs5DisplacedBend")), + nDisplacedTracks_(iConfig.getParameter("nDisplacedTracks")), + dzPVTrk_(iConfig.getParameter("MaxDzTrackPV")), tTopoToken_(esConsumes(edm::ESInputTag("", ""))), trackToken_(consumes>>(iConfig.getParameter("L1TrackInputTag"))), - PVtxToken_(consumes>(iConfig.getParameter("L1PVertexCollection"))) { + PVtxToken_(consumes(iConfig.getParameter("L1PVertexInputTag"))) { zStep_ = 2.0 * trkZMax_ / (zBins_ + 1); // added +1 in denom etaStep_ = 2.0 * trkEtaMax_ / etaBins_; //etaStep is the width of an etabin phiStep_ = 2 * M_PI / phiBins_; ////phiStep is the width of a phibin @@ -140,8 +131,6 @@ L1TrackJetProducer::L1TrackJetProducer(const ParameterSet &iConfig) produces("L1TrackJets"); } -L1TrackJetProducer::~L1TrackJetProducer() {} - void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { unique_ptr L1L1TrackJetProducer(new TkJetCollection); @@ -151,13 +140,14 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { edm::Handle>> TTTrackHandle; iEvent.getByToken(trackToken_, TTTrackHandle); - edm::Handle> PVtx; + edm::Handle PVtx; iEvent.getByToken(PVtxToken_, PVtx); float PVz = (PVtx->at(0)).z0(); L1TrkPtrs_.clear(); tdtrk_.clear(); + // track selection for (unsigned int this_l1track = 0; this_l1track < TTTrackHandle->size(); this_l1track++) { edm::Ptr trkPtr(TTTrackHandle, this_l1track); float trk_pt = trkPtr->momentum().perp(); @@ -165,7 +155,6 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { float trk_chi2dof = trkPtr->chi2Red(); float trk_d0 = trkPtr->d0(); float trk_bendchi2 = trkPtr->stubPtConsistency(); - float trk_z0 = trkPtr->z0(); int trk_nPS = 0; for (int istub = 0; istub < trk_nstubs; istub++) { // loop over the stubs @@ -179,11 +168,22 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { // select tracks if (trk_nPS < trkNPSStubMin_) continue; - if (!trackQualityCuts(trk_nstubs, trk_chi2dof, trk_bendchi2)) + if (!TrackQualitySelection(trk_nstubs, + trk_chi2dof, + trk_bendchi2, + nStubs4PromptBend_, + nStubs5PromptBend_, + nStubs4PromptChi2_, + nStubs5PromptChi2_, + nStubs4DisplacedBend_, + nStubs5DisplacedBend_, + nStubs4DisplacedChi2_, + nStubs5DisplacedChi2_, + displaced_)) continue; - if (std::abs(PVz - trk_z0) > dzPVTrk_ && dzPVTrk_ > 0) + if (std::abs(PVz - trkPtr->z0()) > dzPVTrk_ && dzPVTrk_ > 0) continue; - if (std::abs(trk_z0) > trkZMax_) + if (std::abs(trkPtr->z0()) > trkZMax_) continue; if (std::abs(trkPtr->momentum().eta()) > trkEtaMax_) continue; @@ -198,167 +198,165 @@ void L1TrackJetProducer::produce(Event &iEvent, const EventSetup &iSetup) { tdtrk_.push_back(0); // not displaced track } - if (!L1TrkPtrs_.empty()) { - MaxZBin mzb; - mzb.znum = -1; - mzb.nclust = -1; - mzb.ht = -1; - EtaPhiBin epbins_default[phiBins_][etaBins_]; // create grid of phiBins - - float phi = -1.0 * M_PI; - for (int i = 0; i < phiBins_; ++i) { - float eta = -1.0 * trkEtaMax_; - for (int j = 0; j < etaBins_; ++j) { - epbins_default[i][j].phi = (phi + (phi + phiStep_)) / 2.0; // phimin=phi; phimax= phimin+step - epbins_default[i][j].eta = (eta + (eta + etaStep_)) / 2.0; // phimin=phi; phimax phimin+step - eta += etaStep_; - } // for each etabin - phi += phiStep_; - } // for each phibin (finished creating epbins) - - std::vector zmins, zmaxs; - for (int zbin = 0; zbin < zBins_; zbin++) { - zmins.push_back(-1.0 * trkZMax_ + zStep_ * zbin); - zmaxs.push_back(-1.0 * trkZMax_ + zStep_ * zbin + zStep_ * 2); + // if no tracks pass selection return empty containers + if (L1TrkPtrs_.empty()) { + if (displaced_) + iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJetsExtended"); + else + iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJets"); + return; + } + + MaxZBin mzb; + mzb.znum = -1; + mzb.nclust = -1; + mzb.ht = -1; + + // create 2D grid of eta phi bins + EtaPhiBin epbins_default[phiBins_][etaBins_]; + float phi = -1.0 * M_PI; + for (int i = 0; i < phiBins_; ++i) { + float eta = -1.0 * trkEtaMax_; + for (int j = 0; j < etaBins_; ++j) { + epbins_default[i][j].phi = (phi + (phi + phiStep_)) / 2.0; + epbins_default[i][j].eta = (eta + (eta + etaStep_)) / 2.0; + eta += etaStep_; + } // for each etabin + phi += phiStep_; + } // for each phibin (finished creating bins) + + // create z-bins (might be useful for displaced if we run w/o dz cut) + std::vector zmins, zmaxs; + for (int zbin = 0; zbin < zBins_; zbin++) { + zmins.push_back(-1.0 * trkZMax_ + zStep_ * zbin); + zmaxs.push_back(-1.0 * trkZMax_ + zStep_ * zbin + zStep_ * 2); + } + + // create vectors of clusters + std::vector> L1clusters; + L1clusters.reserve(phiBins_); + std::vector L2clusters; + + // default (empty) grid + for (int i = 0; i < phiBins_; ++i) { + for (int j = 0; j < etaBins_; ++j) { + epbins_default[i][j].pTtot = 0; + epbins_default[i][j].used = false; + epbins_default[i][j].ntracks = 0; + epbins_default[i][j].nxtracks = 0; + epbins_default[i][j].trackidx.clear(); } + } - // create vectors that hold data - std::vector> L1clusters; - L1clusters.reserve(phiBins_); - std::vector L2clusters; - - for (unsigned int zbin = 0; zbin < zmins.size(); ++zbin) { - // initialize matrices - float zmin = zmins[zbin]; - float zmax = zmaxs[zbin]; - EtaPhiBin epbins[phiBins_][etaBins_]; - std::copy(&epbins_default[0][0], &epbins_default[0][0] + phiBins_ * etaBins_, &epbins[0][0]); - - //clear containers - L1clusters.clear(); - L2clusters.clear(); - - // fill grid - for (unsigned int k = 0; k < L1TrkPtrs_.size(); ++k) { - float trkZ = L1TrkPtrs_[k]->z0(); - if (zmax < trkZ) - continue; - if (zmin > trkZ) - continue; - if (zbin == 0 && zmin == trkZ) - continue; - float trkpt = L1TrkPtrs_[k]->momentum().perp(); - float trketa = L1TrkPtrs_[k]->momentum().eta(); - float trkphi = L1TrkPtrs_[k]->momentum().phi(); - for (int i = 0; i < phiBins_; ++i) { - for (int j = 0; j < etaBins_; ++j) { - float eta_min = epbins[i][j].eta - etaStep_ / 2.0; //eta min - float eta_max = epbins[i][j].eta + etaStep_ / 2.0; //eta max - float phi_min = epbins[i][j].phi - phiStep_ / 2.0; //phi min - float phi_max = epbins[i][j].phi + phiStep_ / 2.0; //phi max - - if ((trketa < eta_min) || (trketa > eta_max) || (trkphi < phi_min) || (trkphi > phi_max)) - continue; - - if (trkpt < trkPtMax_) - epbins[i][j].pTtot += trkpt; - else - epbins[i][j].pTtot += trkPtMax_; - epbins[i][j].numtdtrks += tdtrk_[k]; - epbins[i][j].trackidx.push_back(k); - ++epbins[i][j].numtracks; - } // for each phibin - } // for each phibin - } //end loop over tracks - - // first layer clustering - in eta for all phi bins - for (int phibin = 0; phibin < phiBins_; ++phibin) { - L1clusters.push_back(L1_clustering(epbins[phibin], etaBins_, etaStep_)); - } + for (unsigned int zbin = 0; zbin < zmins.size(); ++zbin) { + // initialize grid + float zmin = zmins[zbin]; + float zmax = zmaxs[zbin]; + EtaPhiBin epbins[phiBins_][etaBins_]; + std::copy(&epbins_default[0][0], &epbins_default[0][0] + phiBins_ * etaBins_, &epbins[0][0]); + + //clear cluster containers + L1clusters.clear(); + L2clusters.clear(); + + // fill grid with tracks + for (unsigned int k = 0; k < L1TrkPtrs_.size(); ++k) { + float trkZ = L1TrkPtrs_[k]->z0(); + if (zmax < trkZ) + continue; + if (zmin > trkZ) + continue; + if (zbin == 0 && zmin == trkZ) + continue; + float trkpt = L1TrkPtrs_[k]->momentum().perp(); + float trketa = L1TrkPtrs_[k]->momentum().eta(); + float trkphi = L1TrkPtrs_[k]->momentum().phi(); + for (int i = 0; i < phiBins_; ++i) { + for (int j = 0; j < etaBins_; ++j) { + float eta_min = epbins[i][j].eta - etaStep_ / 2.0; //eta min + float eta_max = epbins[i][j].eta + etaStep_ / 2.0; //eta max + float phi_min = epbins[i][j].phi - phiStep_ / 2.0; //phi min + float phi_max = epbins[i][j].phi + phiStep_ / 2.0; //phi max + + if ((trketa < eta_min) || (trketa > eta_max) || (trkphi < phi_min) || (trkphi > phi_max)) + continue; + + if (trkpt < trkPtMax_) + epbins[i][j].pTtot += trkpt; + else + epbins[i][j].pTtot += trkPtMax_; + epbins[i][j].nxtracks += tdtrk_[k]; + epbins[i][j].trackidx.push_back(k); + ++epbins[i][j].ntracks; + } // for each etabin + } // for each phibin + } //end loop over tracks + + // cluster tracks in eta (first layer) using grid + for (int phibin = 0; phibin < phiBins_; ++phibin) { + L1clusters.push_back(L1_clustering(epbins[phibin], etaBins_, etaStep_)); + } - //second layer clustering in phi for all eta clusters - L2clusters = L2_clustering(L1clusters, phiBins_, phiStep_, etaStep_); - - // sum all cluster pTs in this zbin to find max - float sum_pt = 0; - for (unsigned int k = 0; k < L2clusters.size(); ++k) { - if (L2clusters[k].pTtot > lowpTJetThreshold_ && L2clusters[k].numtracks < lowpTJetMinTrackMultiplicity_) - continue; - if (L2clusters[k].pTtot > highpTJetThreshold_ && L2clusters[k].numtracks < highpTJetMinTrackMultiplicity_) - continue; - if (L2clusters[k].pTtot > minTrkJetpT_) - sum_pt += L2clusters[k].pTtot; - } + // second layer clustering in phi for using eta clusters + L2clusters = L2_clustering(L1clusters, phiBins_, phiStep_, etaStep_); - if (sum_pt < mzb.ht) + // sum all cluster pTs in this zbin to find max + float sum_pt = 0; + for (unsigned int k = 0; k < L2clusters.size(); ++k) { + if (L2clusters[k].pTtot > lowpTJetThreshold_ && L2clusters[k].ntracks < lowpTJetMinTrackMultiplicity_) + continue; + if (L2clusters[k].pTtot > highpTJetThreshold_ && L2clusters[k].ntracks < highpTJetMinTrackMultiplicity_) continue; - // if ht is larger than previous max, this is the new vertex zbin - mzb.ht = sum_pt; - mzb.znum = zbin; - mzb.clusters = L2clusters; - mzb.nclust = L2clusters.size(); - mzb.zbincenter = (zmin + zmax) / 2.0; - } //zbin loop - - vector> L1TrackAssocJet; - for (unsigned int j = 0; j < mzb.clusters.size(); ++j) { - float jetEta = mzb.clusters[j].eta; - float jetPhi = mzb.clusters[j].phi; - float jetPt = mzb.clusters[j].pTtot; - float jetPx = jetPt * cos(jetPhi); - float jetPy = jetPt * sin(jetPhi); - float jetPz = jetPt * sinh(jetEta); - float jetP = jetPt * cosh(jetEta); - int totalDisptrk = mzb.clusters[j].numtdtrks; - bool isDispJet = (totalDisptrk > nDisplacedTracks_ || totalDisptrk == nDisplacedTracks_); - - math::XYZTLorentzVector jetP4(jetPx, jetPy, jetPz, jetP); - L1TrackAssocJet.clear(); - for (unsigned int itrk = 0; itrk < mzb.clusters[j].trackidx.size(); itrk++) - L1TrackAssocJet.push_back(L1TrkPtrs_[mzb.clusters[j].trackidx[itrk]]); - - TkJet trkJet(jetP4, L1TrackAssocJet, mzb.zbincenter, mzb.clusters[j].numtracks, 0, totalDisptrk, 0, isDispJet); - - if (!L1TrackAssocJet.empty()) - L1L1TrackJetProducer->push_back(trkJet); + if (L2clusters[k].pTtot > minTrkJetpT_) + sum_pt += L2clusters[k].pTtot; } - if (displaced_) - iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJetsExtended"); - else - iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJets"); - } -} + if (sum_pt < mzb.ht) + continue; -void L1TrackJetProducer::beginStream(StreamID) {} - -void L1TrackJetProducer::endStream() {} - -bool L1TrackJetProducer::trackQualityCuts(int trk_nstub, float trk_chi2, float trk_bendchi2) { - bool PassQuality = false; - if (!displaced_) { - if (trk_nstub == 4 && trk_bendchi2 < nStubs4PromptBend_ && - trk_chi2 < nStubs4PromptChi2_) // 4 stubs are the lowest track quality and have different cuts - PassQuality = true; - if (trk_nstub > 4 && trk_bendchi2 < nStubs5PromptBend_ && - trk_chi2 < nStubs5PromptChi2_) // above 4 stubs diffent selection imposed (genrally looser) - PassQuality = true; - } else { - if (trk_nstub == 4 && trk_bendchi2 < nStubs4DisplacedBend_ && - trk_chi2 < nStubs4DisplacedChi2_) // 4 stubs are the lowest track quality and have different cuts - PassQuality = true; - if (trk_nstub > 4 && trk_bendchi2 < nStubs5DisplacedBend_ && - trk_chi2 < nStubs5DisplacedChi2_) // above 4 stubs diffent selection imposed (genrally looser) - PassQuality = true; + // if ht is larger than previous max, this is the new vertex zbin + mzb.ht = sum_pt; + mzb.znum = zbin; + mzb.clusters = L2clusters; + mzb.nclust = L2clusters.size(); + mzb.zbincenter = (zmin + zmax) / 2.0; + } //zbin loop + + // output + vector> L1TrackAssocJet; + for (unsigned int j = 0; j < mzb.clusters.size(); ++j) { + float jetEta = mzb.clusters[j].eta; + float jetPhi = mzb.clusters[j].phi; + float jetPt = mzb.clusters[j].pTtot; + float jetPx = jetPt * cos(jetPhi); + float jetPy = jetPt * sin(jetPhi); + float jetPz = jetPt * sinh(jetEta); + float jetP = jetPt * cosh(jetEta); + int totalDisptrk = mzb.clusters[j].nxtracks; + bool isDispJet = (totalDisptrk > nDisplacedTracks_ || totalDisptrk == nDisplacedTracks_); + + math::XYZTLorentzVector jetP4(jetPx, jetPy, jetPz, jetP); + L1TrackAssocJet.clear(); + for (unsigned int itrk = 0; itrk < mzb.clusters[j].trackidx.size(); itrk++) + L1TrackAssocJet.push_back(L1TrkPtrs_[mzb.clusters[j].trackidx[itrk]]); + + TkJet trkJet(jetP4, L1TrackAssocJet, mzb.zbincenter, mzb.clusters[j].ntracks, 0, totalDisptrk, 0, isDispJet); + + L1L1TrackJetProducer->push_back(trkJet); } - return PassQuality; + + std::sort( + L1L1TrackJetProducer->begin(), L1L1TrackJetProducer->end(), [](auto &a, auto &b) { return a.pt() > b.pt(); }); + if (displaced_) + iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJetsExtended"); + else + iEvent.put(std::move(L1L1TrackJetProducer), "L1TrackJets"); } void L1TrackJetProducer::fillDescriptions(ConfigurationDescriptions &descriptions) { - // l1tTrackJets - edm::ParameterSetDescription desc; + ParameterSetDescription desc; desc.add("L1TrackInputTag", edm::InputTag("l1tTTTracksFromTrackletEmulation", "Level1TTTracks")); - desc.add("L1PVertexCollection", edm::InputTag("l1tVertexProducer", "l1vertices")); + desc.add("L1PVertexInputTag", edm::InputTag("l1tVertexFinderEmulator", "l1verticesEmulation")); desc.add("MaxDzTrackPV", 1.0); desc.add("trk_zMax", 15.0); desc.add("trk_ptMax", 200.0); diff --git a/L1Trigger/L1TTrackMatch/python/l1tTrackJetsEmulation_cfi.py b/L1Trigger/L1TTrackMatch/python/l1tTrackJetsEmulation_cfi.py index 97b8651af4bdd..9276870a56a70 100644 --- a/L1Trigger/L1TTrackMatch/python/l1tTrackJetsEmulation_cfi.py +++ b/L1Trigger/L1TTrackMatch/python/l1tTrackJetsEmulation_cfi.py @@ -1,59 +1,46 @@ import FWCore.ParameterSet.Config as cms -l1tTrackJetsEmulation = cms.EDProducer('L1TrackJetEmulationProducer', +l1tTrackJetsEmulation = cms.EDProducer('L1TrackJetEmulatorProducer', L1TrackInputTag= cms.InputTag("l1tGTTInputProducer", "Level1TTTracksConverted"), - VertexInputTag=cms.InputTag("l1tVertexFinderEmulator", "l1verticesEmulation"), - MaxDzTrackPV = cms.double(0.5), + L1PVertexInputTag= cms.InputTag("l1tVertexFinderEmulator","l1verticesEmulation"), + MaxDzTrackPV = cms.double(1.0), trk_zMax = cms.double (15.) , # maximum track z trk_ptMax = cms.double(200.), # maximumum track pT before saturation [GeV] trk_ptMin = cms.double(2.0), # minimum track pt [GeV] trk_etaMax = cms.double(2.4), # maximum track eta - trk_chi2dofMax=cms.double(10.), # maximum track chi2/dof - trk_bendChi2Max=cms.double(2.2), # maximum track bendchi2 + nStubs4PromptChi2=cms.double(10.0), #Prompt track quality flags for loose/tight + nStubs4PromptBend=cms.double(2.2), + nStubs5PromptChi2=cms.double(10.0), + nStubs5PromptBend=cms.double(2.2), trk_nPSStubMin=cms.int32(-1), # minimum PS stubs, -1 means no cut - minTrkJetpT=cms.double(5.), # minimum track pt to be considered for track jet + minTrkJetpT=cms.double(-1.), # minimum track pt to be considered for track jet etaBins=cms.int32(24), phiBins=cms.int32(27), zBins=cms.int32(1), - d0_cutNStubs4=cms.double(0.15), - d0_cutNStubs5=cms.double(0.5), + d0_cutNStubs4=cms.double(-1), + d0_cutNStubs5=cms.double(-1), lowpTJetMinTrackMultiplicity=cms.int32(2), - lowpTJetMinpT=cms.double(50.), + lowpTJetThreshold=cms.double(50.), highpTJetMinTrackMultiplicity=cms.int32(3), - highpTJetMinpT=cms.double(100.), + highpTJetThreshold=cms.double(100.), displaced=cms.bool(False), #Flag for displaced tracks nStubs4DisplacedChi2=cms.double(5.0), #Displaced track quality flags for loose/tight - nStubs4Displacedbend=cms.double(1.7), + nStubs4DisplacedBend=cms.double(1.7), nStubs5DisplacedChi2=cms.double(2.75), - nStubs5Displacedbend=cms.double(3.5), + nStubs5DisplacedBend=cms.double(3.5), nDisplacedTracks=cms.int32(2) #Number of displaced tracks required per jet ) -l1tTrackJetsExtendedEmulation = cms.EDProducer('L1TrackJetEmulationProducer', - L1TrackInputTag= cms.InputTag("l1tGTTInputProducerExtended", "Level1TTTracksExtendedConverted"), - VertexInputTag=cms.InputTag("l1tVertexFinderEmulator", "l1verticesEmulation"), - MaxDzTrackPV = cms.double(4.0), - trk_zMax = cms.double (15.) , # maximum track z - trk_ptMax = cms.double(200.), # maximumum track pT before saturation [GeV] - trk_ptMin = cms.double(3.0), # minimum track pt [GeV] - trk_etaMax = cms.double(2.4), # maximum track eta - trk_chi2dofMax=cms.double(40.), # maximum track chi2/dof - trk_bendChi2Max=cms.double(40.), # maximum track bendchi2 - trk_nPSStubMin=cms.int32(-1), # minimum # PS stubs, -1 means no cut - minTrkJetpT=cms.double(5.), # minimum track pt to be considered for track jet - etaBins=cms.int32(24), - phiBins=cms.int32(27), - zBins=cms.int32(10), - d0_cutNStubs4=cms.double(-1), # -1 excludes nstub=4 from disp tag - d0_cutNStubs5=cms.double(0.22), - lowpTJetMinTrackMultiplicity=cms.int32(2), - lowpTJetMinpT=cms.double(50.), - highpTJetMinTrackMultiplicity=cms.int32(3), - highpTJetMinpT=cms.double(100.), - displaced=cms.bool(True), #Flag for displaced tracks - nStubs4DisplacedChi2=cms.double(3.3), #Disp tracks selection [trkcut excluded diff --git a/L1Trigger/L1TTrackMatch/test/L1TrackObjectNtupleMaker_cfg.py b/L1Trigger/L1TTrackMatch/test/L1TrackObjectNtupleMaker_cfg.py index f8cffe057f9dd..f4f0ea8e8d126 100644 --- a/L1Trigger/L1TTrackMatch/test/L1TrackObjectNtupleMaker_cfg.py +++ b/L1Trigger/L1TTrackMatch/test/L1TrackObjectNtupleMaker_cfg.py @@ -25,8 +25,8 @@ process.load('Configuration.StandardSequences.Services_cff') process.load('Configuration.EventContent.EventContent_cff') process.load('Configuration.StandardSequences.MagneticField_cff') -process.load('Configuration.Geometry.GeometryExtended2026D77Reco_cff') -process.load('Configuration.Geometry.GeometryExtended2026D77_cff') +process.load('Configuration.Geometry.GeometryExtended2026D95Reco_cff') +process.load('Configuration.Geometry.GeometryExtended2026D95_cff') process.load('Configuration.StandardSequences.EndOfProcess_cff') process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') @@ -43,15 +43,7 @@ process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(10)) readFiles = cms.untracked.vstring( - '/store/relval/CMSSW_12_3_0_pre4/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU_123X_mcRun4_realistic_v3_2026D77PU200-v1/2580000/c6df2819-ed05-4b98-8f92-81b7d1b1092e.root', - '/store/relval/CMSSW_12_3_0_pre4/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU_123X_mcRun4_realistic_v3_2026D77PU200-v1/2580000/3f476d95-1ef7-4be6-977b-6bcd1a7c5678.root', - '/store/relval/CMSSW_12_3_0_pre4/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU_123X_mcRun4_realistic_v3_2026D77PU200-v1/2580000/68d651da-4cb7-4bf4-b002-66aecc57a2bc.root', - '/store/relval/CMSSW_12_3_0_pre4/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU_123X_mcRun4_realistic_v3_2026D77PU200-v1/2580000/db0e0ce2-4c5a-4988-9dbd-52066e40b9d2.root', - '/store/relval/CMSSW_12_3_0_pre4/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU_123X_mcRun4_realistic_v3_2026D77PU200-v1/2580000/257a9712-0a96-47b7-897e-f5d980605e46.root', - '/store/relval/CMSSW_12_3_0_pre4/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU_123X_mcRun4_realistic_v3_2026D77PU200-v1/2580000/bee31399-8559-4243-b539-cae1ea897def.root', - '/store/relval/CMSSW_12_3_0_pre4/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU_123X_mcRun4_realistic_v3_2026D77PU200-v1/2580000/24629540-2377-4168-9ae5-518ddd4c43a9.root', - '/store/relval/CMSSW_12_3_0_pre4/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU_123X_mcRun4_realistic_v3_2026D77PU200-v1/2580000/e31ba8f0-332a-4a1a-8bc0-91a12a5fe3db.root', - '/store/relval/CMSSW_12_3_0_pre4/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU_123X_mcRun4_realistic_v3_2026D77PU200-v1/2580000/17902198-4db6-4fcc-9e8c-787991b4db32.root', + '/store/relval/CMSSW_13_0_0/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/130X_mcRun4_realistic_v2_2026D95noPU-v1/00000/16f6615d-f98c-475f-ad33-0e89934b6c7f.root' ) secFiles = cms.untracked.vstring() @@ -67,7 +59,7 @@ useJobReport = cms.untracked.bool(True) ) -process.TFileService = cms.Service("TFileService", fileName = cms.string('GTTObjects_ttbar200PU.root'), closeFileFast = cms.untracked.bool(True)) +process.TFileService = cms.Service("TFileService", fileName = cms.string('GTTObjects_ttbar200PU_v2p2.root'), closeFileFast = cms.untracked.bool(True)) ############################################################ @@ -86,9 +78,8 @@ # DTC emulation -process.load('L1Trigger.TrackerDTC.ProducerES_cff') process.load('L1Trigger.TrackerDTC.ProducerED_cff') -process.dtc = cms.Path(process.TrackerDTCProducer)#*process.TrackerDTCAnalyzer) +process.dtc = cms.Path(process.TrackerDTCProducer) process.load("L1Trigger.TrackFindingTracklet.L1HybridEmulationTracks_cff") process.load("L1Trigger.L1TTrackMatch.l1tTrackSelectionProducer_cfi") @@ -116,8 +107,8 @@ process.l1tTrackFastJets.L1PrimaryVertexTag = cms.InputTag("l1tVertexFinder", "l1vertices") process.l1tTrackFastJetsExtended.L1PrimaryVertexTag = cms.InputTag("l1tVertexFinder", "l1vertices") -process.l1tTrackJets.L1PVertexCollection = cms.InputTag("l1tVertexFinder", "l1vertices") -process.l1tTrackJetsExtended.L1PVertexCollection = cms.InputTag("l1tVertexFinder", "l1vertices") +process.l1tTrackJets.L1PVertexInputTag = cms.InputTag("l1tVertexFinderEmulator","l1verticesEmulation") +process.l1tTrackJetsExtended.L1PVertexInputTag = cms.InputTag("l1tVertexFinderEmulator","l1verticesEmulation") process.l1tTrackerEtMiss.L1VertexInputTag = cms.InputTag("l1tVertexFinder", "l1vertices") process.l1tTrackerHTMiss.L1VertexInputTag = cms.InputTag("l1tVertexFinder", "l1vertices") process.l1tTrackerEtMissExtended.L1VertexInputTag = cms.InputTag("l1tVertexFinder", "l1vertices") @@ -150,7 +141,6 @@ process.pL1GTTInput = cms.Path(process.l1tGTTInputProducerExtended) process.pL1TrackJetsEmu = cms.Path(process.l1tTrackJetsExtendedEmulation) process.pTkMET = cms.Path(process.l1tTrackerEtMissExtended) - #process.pTkMETEmu = cms.Path(process.L1TrackerEmuEtMissExtended) #Doesn't exist process.pTkMHT = cms.Path(process.l1tTrackerHTMissExtended) process.pTkMHTEmulator = cms.Path(process.l1tTrackerEmuHTMissExtended) DISPLACED = 'Displaced'# @@ -235,9 +225,10 @@ process.ntuple = cms.Path(process.L1TrackNtuple) process.out = cms.OutputModule( "PoolOutputModule", - fastCloning = cms.untracked.bool( False ), + outputCommands = process.RAWSIMEventContent.outputCommands, fileName = cms.untracked.string("test.root" ) ) +process.out.outputCommands.append('keep *TTTrack*_*_*_*') process.pOut = cms.EndPath(process.out) @@ -247,5 +238,7 @@ # use this if cluster/stub associators not available # process.schedule = cms.Schedule(process.TTClusterStubTruth,process.TTTracksEmuWithTruth,process.ntuple) -process.schedule = cms.Schedule(process.TTClusterStub, process.TTClusterStubTruth, process.dtc, process.TTTracksEmuWithTruth, process.pL1GTTInput, process.pPV, process.pPVemu, process.pL1TrackSelection, process.pL1TrackJets, process.pL1TrackJetsEmu, -process.pL1TrackFastJets, process.pTkMET, process.pTkMETEmu, process.pTkMHT, process.pTkMHTEmulator, process.ntuple) +process.schedule = cms.Schedule(process.TTClusterStub, process.TTClusterStubTruth, process.dtc, process.TTTracksEmuWithTruth, process.pL1GTTInput, process.pPV, process.pPVemu, process.pL1TrackSelection, process.pL1TrackJets, process.pL1TrackJetsEmu,process.pL1TrackFastJets, process.pTkMET, process.pTkMETEmu, process.pTkMHT, process.pTkMHTEmulator, process.ntuple) + + + diff --git a/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTSAMuonProducer.cc b/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTSAMuonProducer.cc index 5f64b43b10ccc..c6284873729c5 100644 --- a/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTSAMuonProducer.cc +++ b/L1Trigger/Phase2L1GMT/plugins/Phase2L1TGMTSAMuonProducer.cc @@ -127,17 +127,17 @@ void Phase2L1TGMTSAMuonProducer::produce(edm::Event& iEvent, const edm::EventSet // Description: // =========================================================================== SAMuon Phase2L1TGMTSAMuonProducer::Convertl1tMuon(const l1t::Muon& mu, const int bx_) { - ap_uint qual = mu.hwQual(); + qual_sa_t qual = mu.hwQual(); int charge = mu.charge() > 0 ? 0 : 1; - ap_uint pt = round(mu.pt() / LSBpt); - ap_int phi = round(mu.phi() / LSBphi); - ap_int eta = round(mu.eta() / LSBeta); + pt_sa_t pt = round(mu.pt() / LSBpt); + phi_sa_t phi = round(mu.phi() / LSBphi); + eta_sa_t eta = round(mu.eta() / LSBeta); // FIXME: Below are not well defined in phase1 GMT // Using the version from Correlator for now - ap_int z0 = 0; // No tracks info in Phase 1 + z0_sa_t z0 = 0; // No tracks info in Phase 1 // Use 2 bits with LSB = 30cm for BMTF and 25cm for EMTF currently, but subjet to change - ap_int d0 = mu.hwDXY(); + d0_sa_t d0 = mu.hwDXY(); int bstart = 0; wordtype word(0); @@ -148,7 +148,7 @@ SAMuon Phase2L1TGMTSAMuonProducer::Convertl1tMuon(const l1t::Muon& mu, const int bstart = wordconcat(word, bstart, z0, BITSSAZ0); bstart = wordconcat(word, bstart, d0, BITSSAD0); bstart = wordconcat(word, bstart, charge, 1); - bstart = wordconcat(word, bstart, qual, BITSSAQUALITY); + bstart = wordconcat(word, bstart, qual, BITSSAQUAL); SAMuon samuon(mu, charge, pt.to_uint(), eta.to_int(), phi.to_int(), z0.to_int(), d0.to_int(), qual.to_uint()); samuon.setWord(word); diff --git a/L1Trigger/Phase2L1GMT/plugins/TrackMuonMatchAlgorithm.h b/L1Trigger/Phase2L1GMT/plugins/TrackMuonMatchAlgorithm.h index 3df98ee6ed63e..30f3ad6d8a3f5 100644 --- a/L1Trigger/Phase2L1GMT/plugins/TrackMuonMatchAlgorithm.h +++ b/L1Trigger/Phase2L1GMT/plugins/TrackMuonMatchAlgorithm.h @@ -105,7 +105,7 @@ namespace Phase2L1GMT { if (out.size() == maximum) break; l1t::TrackerMuon muon(mu.trkPtr(), mu.charge(), mu.pt(), mu.eta(), mu.phi(), mu.z0(), mu.d0(), mu.quality()); - //muon.setMuonRef(mu.muonRef()); + muon.setMuonRef(mu.muonRef()); for (const auto& stub : mu.stubs()) muon.addStub(stub); out.push_back(muon); @@ -131,9 +131,9 @@ namespace Phase2L1GMT { bstart = 0; bstart = wordconcat(word2, bstart, mu.hwCharge(), 1); - bstart = wordconcat(word2, bstart, mu.hwQual(), BITSGTQUALITY); + bstart = wordconcat(word2, bstart, mu.hwQual(), BITSGTQUAL); bstart = wordconcat(word2, bstart, mu.hwIso(), BITSGTISO); - bstart = wordconcat(word2, bstart, mu.hwBeta(), BITSMUONBETA); + bstart = wordconcat(word2, bstart, mu.hwBeta(), BITSGTBETA); std::array wordout = {{word1, word2}}; mu.setWord(wordout); diff --git a/L1Trigger/Phase2L1GMT/python/gmt_cfi.py b/L1Trigger/Phase2L1GMT/python/gmt_cfi.py index 6e8bcde812de1..e5d24cf7bd955 100644 --- a/L1Trigger/Phase2L1GMT/python/gmt_cfi.py +++ b/L1Trigger/Phase2L1GMT/python/gmt_cfi.py @@ -70,7 +70,7 @@ RelIsoThresholdM = cms.double(0.05), RelIsoThresholdT = cms.double(0.01), verbose = cms.int32(0), - IsodumpForHLS = cms.int32(1), + IsodumpForHLS = cms.int32(0), ), tauto3mu = cms.PSet() diff --git a/L1Trigger/Phase2L1GMT/test/runGMT.py b/L1Trigger/Phase2L1GMT/test/runGMT.py index 4a870ad535b43..26c0aec35ae64 100644 --- a/L1Trigger/Phase2L1GMT/test/runGMT.py +++ b/L1Trigger/Phase2L1GMT/test/runGMT.py @@ -115,30 +115,10 @@ process.dtTriggerPhase2PrimitiveDigis.scenario = 0 process.load("L1Trigger.Phase2L1GMT.gmt_cff") -process.l1tGMTMuons.trackMatching.verbose=1 -process.l1tGMTMuons.verbose=0 -process.l1tGMTMuons.trackConverter.verbose=0 - - - -#process.schedule = cms.Schedule(process.L1TrackTrigger_step,process.pL1TMuonTPS,process.endjob_step,process.e) # Adding MuonTPS - # Path and EndPath definitions process.raw2digi_step = cms.Path(process.RawToDigi) process.L1TrackTrigger_step = cms.Path(process.L1TrackTrigger) -#process.pL1TkPhotonsCrystal = cms.Path(process.L1TkPhotonsCrystal) -#process.pL1TkIsoElectronsCrystal = cms.Path(process.L1TkIsoElectronsCrystal) -#process.pL1TkElectronsLooseCrystal = cms.Path(process.L1TkElectronsLooseCrystal) -#process.pL1TkElectronsLooseHGC = cms.Path(process.L1TkElectronsLooseHGC) -#process.pL1TkElectronsHGC = cms.Path(process.L1TkElectronsHGC) -process.pL1TkMuon = cms.Path(process.L1TkMuons+process.L1TkMuonsTP) -#process.pL1TkElectronsEllipticMatchHGC = cms.Path(process.L1TkElectronsEllipticMatchHGC) -#process.pL1TkElectronsCrystal = cms.Path(process.L1TkElectronsCrystal) -#process.pL1TkPhotonsHGC = cms.Path(process.L1TkPhotonsHGC) -#process.pL1TkIsoElectronsHGC = cms.Path(process.L1TkIsoElectronsHGC) -#process.pL1TkElectronsEllipticMatchCrystal = cms.Path(process.L1TkElectronsEllipticMatchCrystal) -#process.L1simulation_step = cms.Path(process.SimL1Emulator) process.testpath=cms.Path(process.CalibratedDigis*process.dtTriggerPhase2PrimitiveDigis*process.phase2GMT) process.endjob_step = cms.EndPath(process.endOfProcess) process.FEVTDEBUGHLToutput_step = cms.EndPath(process.FEVTDEBUGHLToutput) diff --git a/L1Trigger/Phase2L1ParticleFlow/data/compositeID.json b/L1Trigger/Phase2L1ParticleFlow/data/compositeID.json new file mode 100644 index 0000000000000..53b57c4329f2d --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/data/compositeID.json @@ -0,0 +1 @@ +{"max_depth": 4, "n_trees": 10, "n_classes": 2, "n_features": 11, "trees": [[{"feature": [1, 2, 0, 3, 0, 8, 2, 3, 5, 5, 1, 5, 2, 1, 1, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.140625, 0.264648438, 11.4795256, 3.5, 16.6602802, 8.5, 0.245117188, -2.5, 0.0630929098, 0.585964918, 0.109375, 0.594282091, 0.252929688, 0.171875, 0.171875, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.415820152, 0.556879818, 0.25938049, -0.557010412, -0.523303628, 0.309383601, 0.345235407, -0.177057937, -0.577262163, -0.463166773, -0.295891047, -0.521123528, 0.248503745, -0.140788257, -0.208826631, -0.495121777], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4], "iLeaf": [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]}], [{"feature": [1, 1, 0, 4, 0, 8, 2, 8, 3, 8, 6, 0, 3, 5, 5, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.140625, 0.109375, 12.9026089, -13.5, 15.8960209, 7.5, 0.256835938, 2.5, -2.5, 3.5, 19.5, 8.71144676, 2.5, 0.981220245, 1.01574898, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.358234763, 0.404870719, -0.173584998, 0.459649861, -0.37743035, -0.0318962, 0.344911397, -0.540425003, -0.450546473, -0.387152672, -0.232723728, -0.485472769, 0.149142668, -0.281907678, -0.319123685, -0.469176918], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4], "iLeaf": [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]}], [{"feature": [1, 2, 0, 0, 0, 8, 6, 1, 1, 0, 8, 5, 3, 2, 1, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.140625, 0.256835938, 9.48655319, 6.65358782, 30.484045, 9.5, 16.5, 0.109375, 0.109375, 11.8594818, 9.5, 0.701351345, -2.5, 0.260742188, 0.171875, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.17315872, -0.212719589, 0.398817331, 0.243386999, -0.380143821, -0.169372901, 0.430269539, -0.376156092, -0.378161937, 0.694794357, -0.450359076, -0.180556566, 0.0424509495, -0.275738239, -0.161556199, -0.390238196], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4], "iLeaf": [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]}], [{"feature": [1, 1, 0, 4, 0, 8, 2, 8, 2, 2, 5, 5, 3, 1, 6, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.140625, 0.109375, 7.69071722, -13.5, 14.5519686, 8.5, 0.254882812, 2.5, 0.256835938, 0.219726562, 0.951564014, 0.0752535984, -1.5, 0.171875, 16.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.301308185, 0.304734617, 0.359152883, 0.0685707554, 0.139380425, -0.196517661, 0.26767078, -0.092287831, 0.0782259256, -0.350551367, -0.37530303, -0.175267071, 0.193004847, -0.116770878, -0.22124359, -0.341241091], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4], "iLeaf": [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]}], [{"feature": [1, 3, 1, 3, 5, 3, 3, 5, 4, 10, 5, 3, 5, 3, 3, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.109375, 4.5, 0.171875, -2.5, 0.0630929098, 2.5, 2.5, 0.0503542498, -26.5, 1.5, 0.0889862627, -1.5, 0.0574394017, -2.5, 3.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.528415024, -0.341509908, -0.224257782, 0.332276493, 0.414823771, -0.166062027, -0.161982179, -0.427514702, -0.287400812, 0.144089252, 0.727860332, -0.371313125, -0.353289098, -0.11641188, -0.313981056, -0.364676535], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4], "iLeaf": [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]}], [{"feature": [1, 4, 2, 8, 10, 3, 2, 0, 4, 2, -2, 3, 3, 0, 0, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.109375, -13.5, 0.262695312, 2.5, 6.5, 1.5, 0.290039062, 6.54269695, -39.5, 0.252929688, 0, -2.5, 3.5, 9.02065659, 47.6570053, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.0, -2.0], "children_left": [1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 29, 21, 23, 25, 27, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 30, 22, 24, 26, 28, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.20732594, 0, 0, 0, 0, -0.0745446905, -0.490941644, -0.30292502, 0.351926744, 0.317744046, 0.0299588144, -0.341244638, 0.0837424397, -0.1667739, -0.344421178, -0.299592078, -0.171457857, -0.324574828, 0.0762122124, -1.20732594, -1.20732594], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 11, 11, 12, 12, 13, 13, 14, 14, 10, 10], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4], "iLeaf": [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]}], [{"feature": [1, 0, 2, 5, 2, 0, 0, 0, 8, 10, 3, 8, 5, 5, 6, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.109375, 7.04383469, 0.254882812, 0.0827398598, 0.280273438, 14.939868, 9.47421837, 2.81626415, 9.5, 3.5, 1.5, 3.5, 1.05517507, 0.057422813, 16.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0296402257, 0.372503281, -0.279174328, 0.221156016, 0.313555121, 0.181570679, 0.0492865816, -0.56892544, -0.23718825, -0.00179062958, 0.199230328, -0.261860698, 0.66392386, -0.278425187, -0.0756268799, -0.2803078], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4], "iLeaf": [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]}], [{"feature": [1, 0, 3, 10, 2, 3, 5, 3, 2, 10, 8, 0, 6, 2, 3, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [0.109375, 20.5230389, 2.5, 5.5, 0.280273438, -2.5, 0.0702819526, 4.5, 0.241210938, 3.5, 9.5, 39.9494553, 17.5, 0.291992188, 3.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.183256373, -0.15534921, 0.0191810895, -1.47547853, 0.31037724, 0.204687983, 0.284539312, -0.43629247, -0.310324848, 0.889182568, 0.0579996705, -0.232304126, 0.737877786, -0.190356895, -0.212541878, -0.326160461], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4], "iLeaf": [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]}], [{"feature": [0, 3, 5, 1, 3, 6, 0, -2, 3, 1, 5, 4, 0, 6, 9, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [20.5809975, -2.5, 1.04217935, 0.078125, 1.5, 19.5, 58.5983238, 0, -4.5, 0.234375, 0.0664517879, -10.0, 24.7252579, 11.5, 7.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2.0, -2.0], "children_left": [1, 3, 5, 7, 9, 11, 13, 29, 15, 17, 19, 21, 23, 25, 27, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, 30, 16, 18, 20, 22, 24, 26, 28, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0.363402009, 0, 0, 0, 0, 0, 0, 0, -0.312803566, -0.204151258, 0.0456621125, -0.14882116, 0.20181866, -0.25203824, -1.20155859, 0.277422816, -0.404247284, 0.160260811, 0.580196619, -0.262717366, 0.406997979, -0.171044543, 0.363402009, 0.363402009], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 7, 7], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4], "iLeaf": [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]}], [{"feature": [0, 2, 5, 3, 2, 2, 0, 1, 3, 4, 5, 6, 0, 9, 9, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2], "threshold": [30.4840431, 0.245117188, 1.16153455, -1.5, 0.282226562, 0.280273438, 58.5983238, 0.140625, 2.5, 9.5, 0.948028684, 11.5, 35.0835876, 3.5, 7.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], "children_left": [1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "children_right": [2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1], "value": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0297101662, -0.228048578, 0.0833817348, -0.235787436, -0.063608937, -0.25884518, -0.224756852, -0.356927872, -0.296937495, 0.299284011, 0.524167776, -0.222849742, -0.0959210619, -0.489832073, 0.43333599, -0.154233307], "parent": [-1, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14], "depth": [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4], "iLeaf": [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]}]], "init_predict": [0, 0], "norm": 1} \ No newline at end of file diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h b/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h new file mode 100644 index 0000000000000..bf784f7bdbb5a --- /dev/null +++ b/L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h @@ -0,0 +1,63 @@ +#ifndef CC_INVERSION_H__ +#define CC_INVERSION_H__ + +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 + 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; + } + +} // namespace l1ct +#endif \ No newline at end of file diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/conifer.h b/L1Trigger/Phase2L1ParticleFlow/interface/conifer.h index 61ee6d89b9e85..542806376f964 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/conifer.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/conifer.h @@ -1,8 +1,12 @@ #ifndef CONIFER_CPP_H__ #define CONIFER_CPP_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 */ +#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](auto 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); @@ -152,4 +163,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 b75be595ee9df..560325e46cb0c 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/egamma/pftkegalgo_ref.h @@ -4,6 +4,8 @@ #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" +#include "L1Trigger/Phase2L1ParticleFlow/interface/common/inversion.h" namespace edm { class ParameterSet; @@ -30,6 +32,8 @@ namespace l1ct { std::vector dEtaValues; std::vector dPhiValues; float trkQualityPtMin; // GeV + bool doCompositeTkEle; + unsigned int nCompCandPerCluster; bool writeEgSta; struct IsoParameters { @@ -53,6 +57,18 @@ namespace l1ct { bool doPfIso; EGIsoEleObjEmu::IsoType hwIsoTypeTkEle; EGIsoObjEmu::IsoType hwIsoTypeTkEm; + + struct CompIDParameters { + CompIDParameters(const edm::ParameterSet &); + 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; + int debug = 0; PFTkEGAlgoEmuConfig(const edm::ParameterSet &iConfig); @@ -72,6 +88,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}, @@ -81,6 +99,7 @@ namespace l1ct { bool doPfIso = false, EGIsoEleObjEmu::IsoType hwIsoTypeTkEle = EGIsoEleObjEmu::IsoType::TkIso, EGIsoObjEmu::IsoType hwIsoTypeTkEm = EGIsoObjEmu::IsoType::TkIsoPV, + const CompIDParameters &compIDparams = {-4, 0.214844, "compositeID.json"}, int debug = 0) : nTRACK(nTrack), @@ -99,6 +118,8 @@ namespace l1ct { dEtaValues(dEtaValues), dPhiValues(dPhiValues), trkQualityPtMin(trkQualityPtMin), + doCompositeTkEle(doCompositeTkEle), + nCompCandPerCluster(nCompCandPerCluster), writeEgSta(writeEgSta), tkIsoParams_tkEle(tkIsoParams_tkEle), tkIsoParams_tkEm(tkIsoParams_tkEm), @@ -108,12 +129,13 @@ namespace l1ct { doPfIso(doPfIso), hwIsoTypeTkEle(hwIsoTypeTkEle), hwIsoTypeTkEm(hwIsoTypeTkEm), + compIDparams(compIDparams), debug(debug) {} }; class PFTkEGAlgoEmulator { public: - PFTkEGAlgoEmulator(const PFTkEGAlgoEmuConfig &config) : cfg(config), debug_(cfg.debug) {} + PFTkEGAlgoEmulator(const PFTkEGAlgoEmuConfig &config); virtual ~PFTkEGAlgoEmulator() {} @@ -132,13 +154,33 @@ 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; - 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 +194,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 +208,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 +226,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 +354,7 @@ namespace l1ct { z0_t z0) const; PFTkEGAlgoEmuConfig cfg; + conifer::BDT, false> *composite_bdt_; int debug_; }; } // namespace l1ct diff --git a/L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/hgcalinput_ref.h b/L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/hgcalinput_ref.h index aa99831c369ba..5fd80ba3954d7 100644 --- a/L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/hgcalinput_ref.h +++ b/L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/hgcalinput_ref.h @@ -3,10 +3,18 @@ #include "DataFormats/L1TParticleFlow/interface/layer1_emulator.h" +namespace edm { + class ParameterSet; +} + namespace l1ct { class HgcalClusterDecoderEmulator { + bool slim_; + public: - HgcalClusterDecoderEmulator(){}; + HgcalClusterDecoderEmulator(bool slim = false) : slim_{slim} {}; + HgcalClusterDecoderEmulator(const edm::ParameterSet &pset); + ~HgcalClusterDecoderEmulator(); l1ct::HadCaloObjEmu decode(const ap_uint<256> &in) const; }; 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/L1MHtPFProducer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1MHtPFProducer.cc index 741b7a84e550c..1a970da6e8ad0 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1MHtPFProducer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1MHtPFProducer.cc @@ -65,7 +65,7 @@ void L1MhtPfProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::Even std::vector L1MhtPfProducer::convertEDMToHW(std::vector edmJets) const { std::vector hwJets; std::for_each(edmJets.begin(), edmJets.end(), [&](l1t::PFJet jet) { - l1ct::Jet hwJet = l1ct::Jet::unpack(jet.encodedJet()); + l1ct::Jet hwJet = jet.getHWJetCT(); hwJets.push_back(hwJet); }); return hwJets; @@ -75,17 +75,18 @@ std::vector L1MhtPfProducer::convertHWToEDM(l1ct::Sum hwSums) const std::vector edmSums; reco::Candidate::PolarLorentzVector htVector; - htVector.SetPt(hwSums.hwSumPt.to_double()); + l1gt::Sum gtSum = hwSums.toGT(); + htVector.SetPt(l1gt::Scales::floatPt(gtSum.scalar_pt)); htVector.SetPhi(0); htVector.SetEta(0); reco::Candidate::PolarLorentzVector mhtVector; - mhtVector.SetPt(hwSums.hwPt.to_double()); - mhtVector.SetPhi(l1ct::Scales::floatPhi(hwSums.hwPhi)); + mhtVector.SetPt(l1gt::Scales::floatPt(gtSum.vector_pt)); + mhtVector.SetPhi(l1gt::Scales::floatPhi(gtSum.vector_phi)); mhtVector.SetEta(0); - l1t::EtSum ht(htVector, l1t::EtSum::EtSumType::kTotalHt); - l1t::EtSum mht(mhtVector, l1t::EtSum::EtSumType::kMissingHt); + l1t::EtSum ht(htVector, l1t::EtSum::EtSumType::kTotalHt, gtSum.scalar_pt.bits_to_uint64()); + l1t::EtSum mht(mhtVector, l1t::EtSum::EtSumType::kMissingHt, gtSum.vector_pt.bits_to_uint64(), 0, gtSum.vector_phi); edmSums.push_back(ht); edmSums.push_back(mht); diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1SeedConePFJetProducer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1SeedConePFJetProducer.cc index 4f64649d48103..5cf255385f433 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1SeedConePFJetProducer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1SeedConePFJetProducer.cc @@ -201,7 +201,8 @@ std::vector L1SeedConePFJetProducer::convertHWToEDM( gtJet.v3.pt.V, gtJet.v3.eta.V, gtJet.v3.phi.V); - edmJet.setEncodedJet(jet.toGT().pack()); + edmJet.setEncodedJet(l1t::PFJet::HWEncoding::CT, jet.pack()); + edmJet.setEncodedJet(l1t::PFJet::HWEncoding::GT, jet.toGT().pack()); // get back the references to the constituents std::vector> constituents; std::for_each(jet.constituents.begin(), jet.constituents.end(), [&](auto constituent) { diff --git a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc index c27b0187ef408..d6bd40e646190 100644 --- a/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/plugins/L1TCorrelatorLayer1Producer.cc @@ -21,6 +21,7 @@ #include "L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/tkinput_ref.h" #include "L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/muonGmtToL1ct_ref.h" +#include "L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/hgcalinput_ref.h" #include "L1Trigger/Phase2L1ParticleFlow/interface/regionizer/regionizer_base_ref.h" #include "L1Trigger/Phase2L1ParticleFlow/interface/regionizer/multififo_regionizer_ref.h" #include "L1Trigger/Phase2L1ParticleFlow/interface/regionizer/tdr_regionizer_ref.h" @@ -67,13 +68,13 @@ class L1TCorrelatorLayer1Producer : public edm::stream::EDProducer<> { l1ct::Event event_; std::unique_ptr trackInput_; std::unique_ptr muonInput_; + std::unique_ptr hgcalInput_; std::unique_ptr regionizer_; std::unique_ptr l1pfalgo_; std::unique_ptr l1pualgo_; std::unique_ptr l1tkegalgo_; std::unique_ptr l1tkegsorter_; - bool writeEgSta_; // Region dump const std::string regionDumpName_; bool writeRawHgcalCluster_; @@ -110,11 +111,44 @@ class L1TCorrelatorLayer1Producer : public edm::stream::EDProducer<> { void addRawHgcalCluster(l1ct::DetectorSector> &sec, const l1t::PFCluster &c); + template + void rawHgcalClusterEncode(ap_uint<256> &cwrd, const l1ct::DetectorSector &sec, const l1t::PFCluster &c) const { + cwrd = 0; + 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(); + // 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()); + // 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(); + + 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 + // 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 in the interface) + cwrd(213, 201) = w_srrtot; + cwrd(94, 83) = w_meanz; + cwrd(127, 116) = w_hoe.range(); + } + // fetching outputs std::unique_ptr fetchHadCalo() const; 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, @@ -176,6 +210,7 @@ L1TCorrelatorLayer1Producer::L1TCorrelatorLayer1Producer(const edm::ParameterSet #if 0 // LATER produces("TKVtx"); #endif + produces>("DecodedTK"); for (const auto &tag : iConfig.getParameter>("emClusters")) { emCands_.push_back(consumes(tag)); @@ -200,6 +235,13 @@ L1TCorrelatorLayer1Producer::L1TCorrelatorLayer1Producer(const edm::ParameterSet } else if (muInAlgo != "Ideal") throw cms::Exception("Configuration", "Unsupported muonInputConversionAlgo"); + const std::string &hgcalInAlgo = iConfig.getParameter("hgcalInputConversionAlgo"); + if (hgcalInAlgo == "Emulator") { + hgcalInput_ = std::make_unique( + iConfig.getParameter("hgcalInputConversionParameters")); + } else if (hgcalInAlgo != "Ideal") + throw cms::Exception("Configuration", "Unsupported hgcalInputConversionAlgo"); + const std::string ®algo = iConfig.getParameter("regionizerAlgo"); if (regalgo == "Ideal") { regionizer_ = @@ -372,6 +414,7 @@ void L1TCorrelatorLayer1Producer::produce(edm::Event &iEvent, const edm::EventSe iEvent.put(fetchEmCalo(), "EmCalo"); iEvent.put(fetchHadCalo(), "Calo"); iEvent.put(fetchTracks(), "TK"); + iEvent.put(fetchDecodedTracks(), "DecodedTK"); // Then do the vertexing, and save it out std::vector z0s; @@ -644,7 +687,6 @@ void L1TCorrelatorLayer1Producer::addDecodedTrack(l1ct::DetectorSector &sec, const l1t::PFCluster &c) { l1ct::HadCaloObjEmu calo; - 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 - calo.hwPhi = l1ct::Scales::makePhi(sec.region.localPhi(c.phi())); - calo.hwEmPt = l1ct::Scales::makePtFromFloat(c.emEt()); - calo.hwEmID = c.hwEmID(); + ap_uint<256> word = 0; + rawHgcalClusterEncode(word, sec, c); + if (hgcalInput_) { + calo = hgcalInput_->decode(word); + } else { + 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 + calo.hwPhi = l1ct::Scales::makePhi(sec.region.localPhi(c.phi())); + calo.hwEmPt = l1ct::Scales::makePtFromFloat(c.emEt()); + calo.hwEmID = c.hwEmID(); + calo.hwSrrTot = l1ct::Scales::makeSrrTot(c.sigmaRR()); + calo.hwMeanZ = c.absZBarycenter() < 320. ? l1ct::meanz_t(0) : l1ct::Scales::makeMeanZ(c.absZBarycenter()); + calo.hwHoe = l1ct::Scales::makeHoe(c.hOverE()); + } calo.src = &c; sec.obj.push_back(calo); } void L1TCorrelatorLayer1Producer::addRawHgcalCluster(l1ct::DetectorSector> &sec, const l1t::PFCluster &c) { ap_uint<256> cwrd = 0; - ap_uint<14> w_pt = round(c.pt() / 0.25); - 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(); - - 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; - + rawHgcalClusterEncode(cwrd, sec, c); sec.obj.push_back(cwrd); } void L1TCorrelatorLayer1Producer::addDecodedEmCalo(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 @@ -849,6 +891,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) { @@ -869,6 +942,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) { @@ -879,6 +955,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); } } @@ -1050,6 +1128,7 @@ void L1TCorrelatorLayer1Producer::putEgObjects(edm::Event &iEvent, tkele.setHwQual(egele.hwQual); tkele.setPFIsol(egele.floatRelIso(l1ct::EGIsoEleObjEmu::IsoType::PfIso)); tkele.setEgBinaryWord(egele.pack()); + 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/plugins/PFClusterProducerFromHGC3DClusters.cc b/L1Trigger/Phase2L1ParticleFlow/plugins/PFClusterProducerFromHGC3DClusters.cc index 7d393fe347bd9..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()); @@ -156,6 +157,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 3d1621b655d56..4b267ff987db1 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py @@ -1,6 +1,5 @@ import FWCore.ParameterSet.Config as cms - tkEgAlgoParameters = cms.PSet( nTRACK=cms.uint32(50), # very large numbers for first test nTRACK_EGIN=cms.uint32(50), # very large numbers for first test @@ -50,7 +49,15 @@ doTkIso=cms.bool(True), doPfIso=cms.bool(True), hwIsoTypeTkEle=cms.uint32(0), - hwIsoTypeTkEm=cms.uint32(2) + hwIsoTypeTkEm=cms.uint32(2), + doCompositeTkEle=cms.bool(False), + nCompCandPerCluster=cms.uint32(3), + compositeParametersTkEle=cms.PSet( + # the working points are cuts on BDT output logits log(p/1-p) + loose_wp=cms.double(-0.732422), + tight_wp=cms.double(0.214844), + model=cms.string("L1Trigger/Phase2L1ParticleFlow/data/compositeID.json") + ), ) tkEgSorterParameters = cms.PSet( 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 aa5ebd4147e39..4a0774bdfc3b0 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), @@ -48,6 +49,7 @@ ), muonInputConversionAlgo = cms.string("Emulator"), muonInputConversionParameters = muonInputConversionParameters.clone(), + hgcalInputConversionAlgo = cms.string("Ideal"), regionizerAlgo = cms.string("Ideal"), pfAlgo = cms.string("PFAlgo3"), puAlgo = cms.string("LinearizedPuppi"), @@ -170,6 +172,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), @@ -194,6 +197,10 @@ ), muonInputConversionAlgo = cms.string("Emulator"), muonInputConversionParameters = muonInputConversionParameters.clone(), + hgcalInputConversionAlgo = cms.string("Emulator"), + hgcalInputConversionParameters = cms.PSet( + slim = cms.bool(False) + ), regionizerAlgo = cms.string("Multififo"), regionizerAlgoParameters = cms.PSet( useAlsoVtxCoords = cms.bool(True), @@ -262,7 +269,9 @@ doBremRecovery=True, doEndcapHwQual=True, writeBeforeBremRecovery=False, - writeEGSta=True), + writeEGSta=True, + doCompositeTkEle=True, + trkQualityPtMin=0.), # This should be 10 GeV when doCompositeTkEle=False tkEgSorterAlgo = cms.string("Endcap"), tkEgSorterParameters=tkEgSorterParameters.clone( nObjToSort = 5 @@ -292,6 +301,12 @@ writeRawHgcalCluster = cms.untracked.bool(True) ) +l1tLayer1HGCalElliptic = l1tLayer1HGCal.clone( + tkEgAlgoParameters = l1tLayer1HGCal.tkEgAlgoParameters.clone( + doCompositeTkEle = False, + trkQualityPtMin = 10.) +) + l1tLayer1HGCalExtended = l1tLayer1HGCal.clone(tracks = cms.InputTag('l1tPFTracksFromL1TracksExtended')) l1tLayer1HGCalNoTK = cms.EDProducer("L1TCorrelatorLayer1Producer", @@ -307,6 +322,10 @@ trkPtCut = cms.double(2.0), muonInputConversionAlgo = cms.string("Emulator"), muonInputConversionParameters = muonInputConversionParameters.clone(), + hgcalInputConversionAlgo = cms.string("Emulator"), + hgcalInputConversionParameters = cms.PSet( + slim = cms.bool(True) + ), regionizerAlgo = cms.string("Ideal"), pfAlgo = cms.string("PFAlgoDummy"), puAlgo = cms.string("LinearizedPuppi"), @@ -393,6 +412,7 @@ trkPtCut = cms.double(2.0), muonInputConversionAlgo = cms.string("Ideal"), muonInputConversionParameters = muonInputConversionParameters.clone(), + hgcalInputConversionAlgo = cms.string("Ideal"), regionizerAlgo = cms.string("Ideal"), pfAlgo = cms.string("PFAlgoDummy"), puAlgo = cms.string("LinearizedPuppi"), @@ -530,6 +550,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, @@ -541,10 +605,14 @@ L1TLayer1Task = cms.Task( l1tLayer1Barrel, + l1tLayer1BarrelExtended, l1tLayer1HGCal, + l1tLayer1HGCalExtended, l1tLayer1HGCalNoTK, l1tLayer1HF, l1tLayer1, l1tLayer1Extended, - l1tLayer1EG + l1tLayer1HGCalElliptic, + l1tLayer1EG, + l1tLayer1EGElliptic ) diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer2EG_cff.py b/L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer2EG_cff.py index 517c8e82bf733..e42be9d3b8211 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/python/l1pfJetMet_cff.py b/L1Trigger/Phase2L1ParticleFlow/python/l1pfJetMet_cff.py index a795b5f590952..60d58d2ac23dd 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1pfJetMet_cff.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1pfJetMet_cff.py @@ -25,7 +25,7 @@ phase2_hgcalV11.toModify(_correctedJets, correctorFile = "L1Trigger/Phase2L1ParticleFlow/data/jecs/jecs_20220308.root") from L1Trigger.Phase2L1ParticleFlow.l1tMHTPFProducer_cfi import l1tMHTPFProducer -l1tSCPFL1PuppiCorrectedEmulatorMHT = l1tMHTPFProducer.clone() +l1tSCPFL1PuppiCorrectedEmulatorMHT = l1tMHTPFProducer.clone(jets = 'l1tSCPFL1PuppiCorrectedEmulator') L1TPFJetsTask = cms.Task( l1tLayer2Deregionizer, l1tSCPFL1PF, l1tSCPFL1Puppi, l1tSCPFL1PuppiEmulator, l1tSCPFL1PuppiCorrectedEmulator, l1tSCPFL1PuppiCorrectedEmulatorMHT diff --git a/L1Trigger/Phase2L1ParticleFlow/python/l1tPFTracksFromL1Tracks_cfi.py b/L1Trigger/Phase2L1ParticleFlow/python/l1tPFTracksFromL1Tracks_cfi.py index 35f44e8d32e2e..7ebf7d8236d1c 100644 --- a/L1Trigger/Phase2L1ParticleFlow/python/l1tPFTracksFromL1Tracks_cfi.py +++ b/L1Trigger/Phase2L1ParticleFlow/python/l1tPFTracksFromL1Tracks_cfi.py @@ -26,7 +26,7 @@ ) l1tPFTracksFromL1TracksExtended = l1tPFTracksFromL1Tracks.clone( - L1TrackTag = cms.InputTag("TTTracksFromExtendedTrackletEmulation", "Level1TTTracks"), + L1TrackTag = cms.InputTag("l1tTTTracksFromExtendedTrackletEmulation", "Level1TTTracks"), nParam = 5, qualityBits = cms.vstring( "momentum.perp > 2 && getStubRefs.size >= 4 && chi2Red < 15 && POCA.x < 1.0 && POCA.x > -1.0 && POCA.y < 1.0 && POCA.y > -1.0", diff --git a/L1Trigger/Phase2L1ParticleFlow/src/CaloClusterer.cc b/L1Trigger/Phase2L1ParticleFlow/src/CaloClusterer.cc index 269bcc7de53c0..c42d8614a0128 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/CaloClusterer.cc +++ b/L1Trigger/Phase2L1ParticleFlow/src/CaloClusterer.cc @@ -589,6 +589,8 @@ void l1tpf_calo::SimpleCaloLinker::run() { float wecal = cluster.ecal_et / cluster.et, whcal = 1.0 - wecal; cluster.eta = ecal.eta * wecal + hcal.eta * whcal; cluster.phi = ecal.phi * wecal + hcal.phi * whcal; + cluster.ecal_eta = cluster.eta; + cluster.ecal_phi = cluster.phi; // wrap around phi cluster.phi = reco::reduceRange(cluster.phi); cluster.constituents.emplace_back(-i - 1, 1); @@ -617,6 +619,8 @@ void l1tpf_calo::SimpleCaloLinker::run() { cluster.et = myet + etot; cluster.eta = hcal.eta + avg_eta / cluster.et; cluster.phi = hcal.phi + avg_phi / cluster.et; + cluster.ecal_eta = cluster.eta; + cluster.ecal_phi = cluster.phi; // wrap around phi cluster.phi = reco::reduceRange(cluster.phi); } @@ -676,6 +680,8 @@ void l1tpf_calo::FlatCaloLinker::run() { dst.et = src.et; dst.eta = src.eta; dst.phi = src.phi; + dst.ecal_eta = src.eta; + dst.ecal_phi = src.phi; dst.ecal_et = 0; dst.hcal_et = 0; for (const auto &pair : src.constituents) { 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 89ee83ee6b481..da9baff109831 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/egamma/pftkegalgo_ref.cpp @@ -7,6 +7,7 @@ #include #include #include +#include using namespace l1ct; @@ -30,6 +31,8 @@ l1ct::PFTkEGAlgoEmuConfig::PFTkEGAlgoEmuConfig(const edm::ParameterSet &pset) dEtaValues(pset.getParameter>("dEtaValues")), dPhiValues(pset.getParameter>("dPhiValues")), trkQualityPtMin(pset.getParameter("trkQualityPtMin")), + doCompositeTkEle(pset.getParameter("doCompositeTkEle")), + nCompCandPerCluster(pset.getParameter("nCompCandPerCluster")), writeEgSta(pset.getParameter("writeEGSta")), tkIsoParams_tkEle(pset.getParameter("tkIsoParametersTkEle")), tkIsoParams_tkEm(pset.getParameter("tkIsoParametersTkEm")), @@ -39,6 +42,7 @@ l1ct::PFTkEGAlgoEmuConfig::PFTkEGAlgoEmuConfig(const edm::ParameterSet &pset) doPfIso(pset.getParameter("doPfIso")), hwIsoTypeTkEle(static_cast(pset.getParameter("hwIsoTypeTkEle"))), hwIsoTypeTkEm(static_cast(pset.getParameter("hwIsoTypeTkEm"))), + compIDparams(pset.getParameter("compositeParametersTkEle")), debug(pset.getUntrackedParameter("debug", 0)) {} l1ct::PFTkEGAlgoEmuConfig::IsoParameters::IsoParameters(const edm::ParameterSet &pset) @@ -47,8 +51,25 @@ l1ct::PFTkEGAlgoEmuConfig::IsoParameters::IsoParameters(const edm::ParameterSet pset.getParameter("dRMin"), pset.getParameter("dRMax")) {} +l1ct::PFTkEGAlgoEmuConfig::CompIDParameters::CompIDParameters(const edm::ParameterSet &pset) + : CompIDParameters(pset.getParameter("loose_wp"), + pset.getParameter("tight_wp"), + pset.getParameter("model")) {} + #endif +PFTkEGAlgoEmulator::PFTkEGAlgoEmulator(const PFTkEGAlgoEmuConfig &config) + : cfg(config), composite_bdt_(nullptr), debug_(cfg.debug) { + if (cfg.doCompositeTkEle) { +#ifdef CMSSW_GIT_HASH + auto resolvedFileName = edm::FileInPath(cfg.compIDparams.conifer_model).fullPath(); +#else + auto resolvedFileName = cfg.compIDparams.conifer_model; +#endif + composite_bdt_ = new conifer::BDT(resolvedFileName); + } +} + void PFTkEGAlgoEmulator::toFirmware(const PFInputRegion &in, PFRegion ®ion, EmCaloObj emcalo[/*nCALO*/], @@ -109,10 +130,10 @@ 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]; @@ -156,6 +177,90 @@ 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 { + 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 dR2 = (d_phi * d_phi) + (d_eta * d_eta); + + 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; + cand.track_idx = itk; + cand.dpt = std::abs(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(), cfg.nCompCandPerCluster); + if (nCandPerCluster == 0) + continue; + + float maxScore = -999; + int ibest = -1; + for (unsigned int icand = 0; icand < nCandPerCluster; icand++) { + auto &cand = candidates[icand]; + const std::vector &emcalo_sel = emcalo; + float score = compute_composite_score(cand, emcalo_sel, track, cfg.compIDparams); + if ((score > cfg.compIDparams.bdtScore_loose_wp) && (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]; + + // 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_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; + 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 = {tkpt, hoe, srrtot, deta, dphi, dpt, meanz, nstubs, chi2rphi, chi2rz, chi2bend}; + std::vector bdt_score = composite_bdt_->decision_function(inputs); + + return bdt_score[0]; +} + void PFTkEGAlgoEmulator::sel_emCalo(unsigned int nmax_sel, const std::vector &emcalo, std::vector &emcalo_sel) const { @@ -181,7 +286,6 @@ void PFTkEGAlgoEmulator::run(const PFInputRegion &in, OutputRegion &out) const { << std::endl; } } - // FIXME: can be removed in the endcap since now running with the "interceptor". // Might still be needed in barrel // filter and select first N elements of input clusters @@ -193,12 +297,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()); @@ -215,6 +325,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 { @@ -230,6 +341,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) { @@ -242,7 +354,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) @@ -265,7 +377,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); } } @@ -319,7 +431,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; @@ -327,9 +440,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; @@ -338,6 +456,7 @@ EGIsoEleObjEmu &PFTkEGAlgoEmulator::addEGIsoEleToPF(std::vector egiso.hwCharge = track.hwCharge; egiso.srcCluster = calo.src; egiso.srcTrack = track.src; + egiso.idScore = bdtScore; egobjs.push_back(egiso); if (debug_ > 2) @@ -356,6 +475,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()) { @@ -365,7 +485,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/l1-converters/hgcalinputt_ref.cpp b/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp index d1777586ca753..e4713944edf02 100644 --- a/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp +++ b/L1Trigger/Phase2L1ParticleFlow/src/l1-converters/hgcalinputt_ref.cpp @@ -1,5 +1,12 @@ #include "L1Trigger/Phase2L1ParticleFlow/interface/l1-converters/hgcalinput_ref.h" +#ifdef CMSSW_GIT_HASH +#include "FWCore/ParameterSet/interface/ParameterSet.h" +l1ct::HgcalClusterDecoderEmulator::HgcalClusterDecoderEmulator(const edm::ParameterSet &pset) + : slim_(pset.getParameter("slim")) {} + +#endif + l1ct::HgcalClusterDecoderEmulator::~HgcalClusterDecoderEmulator() {} l1ct::HadCaloObjEmu l1ct::HgcalClusterDecoderEmulator::decode(const ap_uint<256> &in) const { @@ -8,6 +15,10 @@ 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(); @@ -16,6 +27,11 @@ l1ct::HadCaloObjEmu l1ct::HgcalClusterDecoderEmulator::decode(const ap_uint<256> out.hwPhi = w_phi; // relative to the region center, at calo out.hwEmPt = w_empt * l1ct::pt_t(l1ct::Scales::INTPT_LSB); out.hwEmID = w_qual; - + if (!slim_) { + 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); + } 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_l1ct_binaryFiles_cfg.py b/L1Trigger/Phase2L1ParticleFlow/test/make_l1ct_binaryFiles_cfg.py index 1ad6414ee65bb..29f94aed7b5ba 100644 --- a/L1Trigger/Phase2L1ParticleFlow/test/make_l1ct_binaryFiles_cfg.py +++ b/L1Trigger/Phase2L1ParticleFlow/test/make_l1ct_binaryFiles_cfg.py @@ -119,6 +119,7 @@ process.l1tLayer1BarrelSerenity + process.l1tLayer1Barrel9 + process.l1tLayer1HGCal + + process.l1tLayer1HGCalElliptic + process.l1tLayer1HGCalNoTK + process.l1tLayer1HF + process.l1tLayer1 + @@ -146,7 +147,7 @@ process.l1tLayer2SeedConeJetWriter.maxLinesPerFile = eventsPerFile_*54 if not args.dumpFilesOFF: - for det in "Barrel", "BarrelSerenity", "Barrel9", "HGCal", "HGCalNoTK", "HF": + for det in "Barrel", "BarrelSerenity", "Barrel9", "HGCal", "HGCalElliptic", "HGCalNoTK", "HF": l1pf = getattr(process, 'l1tLayer1'+det) l1pf.dumpFileName = cms.untracked.string("TTbar_PU200_"+det+".dump") diff --git a/L1Trigger/TrackFindingTracklet/test/L1TrackNtupleMaker_cfg.py b/L1Trigger/TrackFindingTracklet/test/L1TrackNtupleMaker_cfg.py index 17a89466a4c2f..eeb0ea60cd8ac 100644 --- a/L1Trigger/TrackFindingTracklet/test/L1TrackNtupleMaker_cfg.py +++ b/L1Trigger/TrackFindingTracklet/test/L1TrackNtupleMaker_cfg.py @@ -11,9 +11,14 @@ # edit options here ############################################################ -GEOMETRY = "D49" -# Set L1 tracking algorithm: -# 'HYBRID' (baseline, 4par fit) or 'HYBRID_DISPLACED' (extended, 5par fit). +# D76 used for old CMSSW_11_3 MC datasets. D88 used for CMSSW_12_6 datasets. +#GEOMETRY = "D76" +GEOMETRY = "D88" + +# Set L1 tracking algorithm: +# 'HYBRID' (baseline, 4par fit) or 'HYBRID_DISPLACED' (extended, 5par fit). +# 'HYBRID_NEWKF' (baseline, 4par fit, with bit-accurate KF emulation), +# 'HYBRID_REDUCED' to use the "Summer Chain" configuration with reduced inputs. # (Or legacy algos 'TMTT' or 'TRACKLET'). L1TRKALGO = 'HYBRID' @@ -31,14 +36,12 @@ process.MessageLogger.L1track = dict(limit = -1) process.MessageLogger.Tracklet = dict(limit = -1) -if GEOMETRY == "D49": - print "using geometry " + GEOMETRY + " (tilted)" - process.load('Configuration.Geometry.GeometryExtended2026D49Reco_cff') - process.load('Configuration.Geometry.GeometryExtended2026D49_cff') -elif GEOMETRY == "D76": - print "using geometry " + GEOMETRY + " (tilted)" - process.load('Configuration.Geometry.GeometryExtended2026D76Reco_cff') - process.load('Configuration.Geometry.GeometryExtended2026D76_cff') +if GEOMETRY == "D76" or GEOMETRY == "D88": + # Use D88 for both, as both correspond to identical CMS Tracker design, and D76 + # unavailable in CMSSW_12_6_0. + print("using geometry " + GEOMETRY + " (tilted)") + process.load('Configuration.Geometry.GeometryExtended2026D88Reco_cff') + process.load('Configuration.Geometry.GeometryExtended2026D88_cff') else: print "this is not a valid geometry!!!" @@ -61,13 +64,13 @@ #from MCsamples.Scripts.getCMSdata_cfi import * #from MCsamples.Scripts.getCMSlocaldata_cfi import * -if GEOMETRY == "D49": +if GEOMETRY == "D76": # Read data from card files (defines getCMSdataFromCards()): #from MCsamples.RelVal_1120.PU200_TTbar_14TeV_cfi import * #inputMC = getCMSdataFromCards() # Or read .root files from directory on local computer: - #dirName = "$myDir/whatever/" + #dirName = "$scratchmc/MCsamples1130_D76/RelVal/TTbar/PU200/" #inputMC=getCMSlocaldata(dirName) # Or read specified dataset (accesses CMS DB, so use this method only occasionally): @@ -79,11 +82,30 @@ elif GEOMETRY == "D76": inputMC = ["/store/relval/CMSSW_11_3_0_pre6/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU_113X_mcRun4_realistic_v6_2026D76PU200-v1/00000/00026541-6200-4eed-b6f8-d3a1fd720e9c.root"] + +elif GEOMETRY == "D88": + + # Read specified .root file: + inputMC = ["/store/relval/CMSSW_12_6_0_pre4/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU_125X_mcRun4_realistic_v2_2026D88PU200-v1/2590000/00b3d04b-4c7b-4506-8d82-9538fb21ee19.root"] + else: print "this is not a valid geometry!!!" process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring(*inputMC)) +if GEOMETRY == "D76": + # If reading old MC dataset, drop incompatible EDProducts. + process.source.dropDescendantsOfDroppedBranches = cms.untracked.bool(False) + process.source.inputCommands = cms.untracked.vstring() + process.source.inputCommands.append('keep *_*_*Level1TTTracks*_*') + process.source.inputCommands.append('keep *_*_*StubAccepted*_*') + process.source.inputCommands.append('keep *_*_*ClusterAccepted*_*') + process.source.inputCommands.append('keep *_*_*MergedTrackTruth*_*') + process.source.inputCommands.append('keep *_genParticles_*_*') + +# Use skipEvents to select particular single events for test vectors +#process.source.skipEvents = cms.untracked.uint32(11) + process.TFileService = cms.Service("TFileService", fileName = cms.string('TTbar_PU200_'+GEOMETRY+'.root'), closeFileFast = cms.untracked.bool(True)) process.Timing = cms.Service("Timing", summaryOnly = cms.untracked.bool(True)) @@ -93,6 +115,8 @@ ############################################################ process.load('L1Trigger.TrackTrigger.TrackTrigger_cff') +process.load('L1Trigger.TrackerTFP.ProducerES_cff') +process.load('L1Trigger.TrackerTFP.ProducerLayerEncoding_cff') # remake stubs? #from L1Trigger.TrackTrigger.TTStubAlgorithmRegister_cfi import * @@ -142,11 +166,32 @@ L1TRK_NAME = "l1tTTTracksFromExtendedTrackletEmulation" L1TRK_LABEL = "Level1TTTracks" L1TRUTH_NAME = "TTTrackAssociatorFromPixelDigisExtended" - -# LEGACY ALGORITHM (EXPERTS ONLY): TRACKLET + +# HYBRID_NEWKF: prompt tracking or reduced +elif (L1TRKALGO == 'HYBRID_NEWKF' or L1TRKALGO == 'HYBRID_REDUCED'): + process.load( 'L1Trigger.TrackFindingTracklet.Producer_cff' ) + NHELIXPAR = 4 + L1TRK_NAME = process.TrackFindingTrackletProducer_params.LabelTT.value() + L1TRK_LABEL = process.TrackFindingTrackletProducer_params.BranchAcceptedTracks.value() + L1TRUTH_NAME = "TTTrackAssociatorFromPixelDigis" + process.TTTrackAssociatorFromPixelDigis.TTTracks = cms.VInputTag( cms.InputTag(L1TRK_NAME, L1TRK_LABEL) ) + process.HybridNewKF = cms.Sequence(process.L1THybridTracks + process.TrackFindingTrackletProducerTBout + process.TrackFindingTrackletProducerKFin + process.TrackFindingTrackletProducerKF + process.TrackFindingTrackletProducerTT + process.TrackFindingTrackletProducerAS + process.TrackFindingTrackletProducerKFout) + process.TTTracksEmulation = cms.Path(process.HybridNewKF) + #process.TTTracksEmulationWithTruth = cms.Path(process.HybridNewKF + process.TrackTriggerAssociatorTracks) + # Optionally include code producing performance plots & end-of-job summary. + process.load( 'SimTracker.TrackTriggerAssociation.StubAssociator_cff' ) + process.load( 'L1Trigger.TrackFindingTracklet.Analyzer_cff' ) + process.TTTracksEmulationWithTruth = cms.Path(process.HybridNewKF + process.TrackTriggerAssociatorTracks + process.StubAssociator + process.TrackFindingTrackletAnalyzerTracklet + process.TrackFindingTrackletAnalyzerTBout + process.TrackFindingTrackletAnalyzerKFin + process.TrackFindingTrackletAnalyzerKF + process.TrackFindingTrackletAnalyzerKFout) + from L1Trigger.TrackFindingTracklet.Customize_cff import * + if (L1TRKALGO == 'HYBRID_NEWKF'): + fwConfig( process ) + if (L1TRKALGO == 'HYBRID_REDUCED'): + reducedConfig( process ) + +# LEGACY ALGORITHM (EXPERTS ONLY): TRACKLET elif (L1TRKALGO == 'TRACKLET'): - print "\n WARNING: This is not the baseline algorithm! Prefer HYBRID or HYBRID_DISPLACED!" - print "\n To run the Tracklet-only algorithm, ensure you have commented out 'CXXFLAGS=-DUSEHYBRID' in BuildFile.xml & recompiled! \n" + print("\n WARNING: This is not the baseline algorithm! Prefer HYBRID or HYBRID_DISPLACED!") + print("\n To run the Tracklet-only algorithm, ensure you have commented out 'CXXFLAGS=-DUSEHYBRID' in BuildFile.xml & recompiled! \n") process.TTTracksEmulation = cms.Path(process.L1THybridTracks) process.TTTracksEmulationWithTruth = cms.Path(process.L1THybridTracksWithAssociators) NHELIXPAR = 4 @@ -248,4 +293,6 @@ process.pd = cms.EndPath(process.writeDataset) process.schedule.append(process.pd) - + + +