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 0961b61689418..a597e588a84bc 100644 --- a/DataFormats/L1Trigger/src/classes_def.xml +++ b/DataFormats/L1Trigger/src/classes_def.xml @@ -305,7 +305,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/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 1294e4fd1b563..0000000000000 --- a/L1Trigger/L1TTrackMatch/interface/L1Clustering.h +++ /dev/null @@ -1,251 +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--; - } // 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 b42a7dd599c7b..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 = 0; //Numbered from 0 to nzbins (16, 32, or 64) in order - int nclust = 0; //number of clusters in this bin - z0_intern zbincenter = 0; - TrackJetEmulationEtaPhiBin *clusters = nullptr; //list of all the clusters in this bin - pt_intern ht = 0; //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 f4cbc302956af..0000000000000 --- a/L1Trigger/L1TTrackMatch/plugins/L1TrackJetEmulationProducer.cc +++ /dev/null @@ -1,680 +0,0 @@ -/* Software to emulate the hardware 2-layer jet-finding algorithm (fixed-point). *Layers 1 and 2* - * - * Created based on Rishi Patel's L1TrackJetProducer.cc file. - * Authors: Samuel Edwin Leigh, Tyler Wu - * Rutgers, the State University of New Jersey - * Modifications: Georgios Karathanasis - * georgios.karathanasis@cern.ch, CU Boulder - * - * Last update: 19-9-2022 (by GK) - */ - -//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; - //values are initialized by 0 as struct assigns default values for members - TrackJetEmulationMaxZBin all_zBins[nz]; - - 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/TrackFindingTracklet/test/L1TrackNtupleMaker_cfg.py b/L1Trigger/TrackFindingTracklet/test/L1TrackNtupleMaker_cfg.py index 7bf733e0277c3..0d27aea2ee548 100644 --- a/L1Trigger/TrackFindingTracklet/test/L1TrackNtupleMaker_cfg.py +++ b/L1Trigger/TrackFindingTracklet/test/L1TrackNtupleMaker_cfg.py @@ -11,7 +11,10 @@ # edit options here ############################################################ -GEOMETRY = "D76" +# 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), @@ -34,14 +37,12 @@ process.MessageLogger.Tracklet = dict(limit = -1) process.MessageLogger.TrackTriggerHPH = 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": +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.GeometryExtended2026D76Reco_cff') - process.load('Configuration.Geometry.GeometryExtended2026D76_cff') + process.load('Configuration.Geometry.GeometryExtended2026D88Reco_cff') + process.load('Configuration.Geometry.GeometryExtended2026D88_cff') else: print("this is not a valid geometry!!!") @@ -64,16 +65,13 @@ #from MCsamples.Scripts.getCMSdata_cfi import * #from MCsamples.Scripts.getCMSlocaldata_cfi import * -if GEOMETRY == "D49": - inputMC = ["/store/relval/CMSSW_11_3_0_pre3/RelValTTbar_14TeV/GEN-SIM-DIGI-RAW/PU_113X_mcRun4_realistic_v3_2026D49PU200_rsb-v1/00000/00260a30-734a-4a3a-a4b0-f836ce5502c6.root"] - -elif GEOMETRY == "D76": +if GEOMETRY == "D76": # Read data from card files (defines getCMSdataFromCards()): #from MCsamples.RelVal_1130_D76.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): @@ -83,13 +81,29 @@ # Or read specified .root file: 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 = cms.Source("PoolSource", fileNames = cms.untracked.vstring(*inputMC), skipEvents = cms.untracked.uint32(11)) +#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)) @@ -100,6 +114,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 * @@ -157,7 +173,7 @@ 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.L1HybridTracks + process.TrackFindingTrackletProducerTBout + process.TrackFindingTrackletProducerKFin + process.TrackFindingTrackletProducerKF + process.TrackFindingTrackletProducerTT + process.TrackFindingTrackletProducerAS + process.TrackFindingTrackletProducerKFout) + 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. @@ -172,8 +188,8 @@ # 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 @@ -274,3 +290,6 @@ process.pd = cms.EndPath(process.writeDataset) process.schedule.append(process.pd) + + +