File indexing completed on 2023-05-05 02:47:48
0001 #ifndef L1Trigger_L1TTrackMatch_L1TrackJetClustering_HH
0002 #define L1Trigger_L1TTrackMatch_L1TrackJetClustering_HH
0003 #include <iostream>
0004 #include <fstream>
0005 #include <cmath>
0006 #include <cstdlib>
0007 #include <string>
0008 #include <cstdlib>
0009 #include "DataFormats/L1Trigger/interface/TkJetWord.h"
0010 #include "DataFormats/L1TrackTrigger/interface/TTTrack_TrackWord.h"
0011 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 namespace l1ttrackjet {
0015
0016 const unsigned int PT_INTPART_BITS{9};
0017 const unsigned int ETA_INTPART_BITS{3};
0018 const unsigned int kExtraGlobalPhiBit{4};
0019
0020 typedef ap_ufixed<TTTrack_TrackWord::TrackBitWidths::kRinvSize - 1, PT_INTPART_BITS, AP_TRN, AP_SAT> pt_intern;
0021 typedef ap_fixed<TTTrack_TrackWord::TrackBitWidths::kTanlSize, ETA_INTPART_BITS, AP_TRN, AP_SAT> glbeta_intern;
0022 typedef ap_int<TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit> glbphi_intern;
0023 typedef ap_int<TTTrack_TrackWord::TrackBitWidths::kZ0Size> z0_intern;
0024 typedef ap_uint<TTTrack_TrackWord::TrackBitWidths::kD0Size> d0_intern;
0025
0026 inline const unsigned int DoubleToBit(double value, unsigned int maxBits, double step) {
0027 unsigned int digitized_value = std::floor(std::abs(value) / step);
0028 unsigned int digitized_maximum = (1 << (maxBits - 1)) - 1;
0029 if (digitized_value > digitized_maximum)
0030 digitized_value = digitized_maximum;
0031 if (value < 0)
0032 digitized_value = (1 << maxBits) - digitized_value;
0033 return digitized_value;
0034 }
0035 inline const double BitToDouble(unsigned int bits, unsigned int maxBits, double step) {
0036 int isign = 1;
0037 unsigned int digitized_maximum = (1 << maxBits) - 1;
0038 if (bits & (1 << (maxBits - 1))) {
0039 isign = -1;
0040 bits = (1 << (maxBits + 1)) - bits;
0041 }
0042 return (double(bits & digitized_maximum) + 0.5) * step * isign;
0043 }
0044
0045
0046 struct EtaPhiBin {
0047 float pTtot;
0048 int ntracks;
0049 int nxtracks;
0050 bool used;
0051 float phi;
0052 float eta;
0053 std::vector<unsigned int> trackidx;
0054 };
0055
0056 struct MaxZBin {
0057 int znum;
0058 int nclust;
0059 float zbincenter;
0060 std::vector<EtaPhiBin> clusters;
0061 float ht;
0062 };
0063
0064
0065 struct TrackJetEmulationEtaPhiBin {
0066 pt_intern pTtot;
0067 l1t::TkJetWord::nt_t ntracks;
0068 l1t::TkJetWord::nx_t nxtracks;
0069 bool used;
0070 glbphi_intern phi;
0071 glbeta_intern eta;
0072 std::vector<unsigned int> trackidx;
0073 };
0074
0075
0076 struct TrackJetEmulationMaxZBin {
0077 int znum;
0078 int nclust;
0079 z0_intern zbincenter;
0080 std::vector<TrackJetEmulationEtaPhiBin> clusters;
0081 pt_intern ht;
0082 };
0083
0084
0085 inline bool TrackQualitySelection(int trk_nstub,
0086 double trk_chi2,
0087 double trk_bendchi2,
0088 double nStubs4PromptBend_,
0089 double nStubs5PromptBend_,
0090 double nStubs4PromptChi2_,
0091 double nStubs5PromptChi2_,
0092 double nStubs4DisplacedBend_,
0093 double nStubs5DisplacedBend_,
0094 double nStubs4DisplacedChi2_,
0095 double nStubs5DisplacedChi2_,
0096 bool displaced_) {
0097 bool PassQuality = false;
0098 if (!displaced_) {
0099 if (trk_nstub == 4 && trk_bendchi2 < nStubs4PromptBend_ &&
0100 trk_chi2 < nStubs4PromptChi2_)
0101 PassQuality = true;
0102 if (trk_nstub > 4 && trk_bendchi2 < nStubs5PromptBend_ &&
0103 trk_chi2 < nStubs5PromptChi2_)
0104 PassQuality = true;
0105 } else {
0106 if (trk_nstub == 4 && trk_bendchi2 < nStubs4DisplacedBend_ &&
0107 trk_chi2 < nStubs4DisplacedChi2_)
0108 PassQuality = true;
0109 if (trk_nstub > 4 && trk_bendchi2 < nStubs5DisplacedBend_ &&
0110 trk_chi2 < nStubs5DisplacedChi2_)
0111 PassQuality = true;
0112 }
0113 return PassQuality;
0114 }
0115
0116
0117 template <typename T, typename Pt, typename Eta, typename Phi>
0118 inline std::vector<T> L1_clustering(T *phislice, int etaBins_, Eta etaStep_) {
0119 std::vector<T> clusters;
0120
0121 int nclust = 0;
0122
0123
0124 for (int etabin = 0; etabin < etaBins_; ++etabin) {
0125 Pt my_pt = 0;
0126 Pt previousbin_pt = 0;
0127 Pt nextbin_pt = 0;
0128 Pt nextbin2_pt = 0;
0129
0130
0131 if (phislice[etabin].used)
0132 continue;
0133
0134 my_pt = phislice[etabin].pTtot;
0135 if (my_pt == 0)
0136 continue;
0137
0138
0139 if (etabin > 0 && !phislice[etabin - 1].used)
0140 previousbin_pt = phislice[etabin - 1].pTtot;
0141
0142
0143 if (etabin < etaBins_ - 1 && !phislice[etabin + 1].used) {
0144 nextbin_pt = phislice[etabin + 1].pTtot;
0145 if (etabin < etaBins_ - 2 && !phislice[etabin + 2].used) {
0146 nextbin2_pt = phislice[etabin + 2].pTtot;
0147 }
0148 }
0149
0150 if (my_pt < previousbin_pt || my_pt <= nextbin_pt) {
0151
0152 if (previousbin_pt > 0) {
0153 clusters.push_back(phislice[etabin - 1]);
0154 phislice[etabin - 1].used = true;
0155 nclust++;
0156 }
0157 continue;
0158 }
0159
0160 clusters.push_back(phislice[etabin]);
0161 phislice[etabin].used = true;
0162 if (previousbin_pt > 0) {
0163 clusters[nclust].pTtot += previousbin_pt;
0164 clusters[nclust].ntracks += phislice[etabin - 1].ntracks;
0165 clusters[nclust].nxtracks += phislice[etabin - 1].nxtracks;
0166 for (unsigned int itrk = 0; itrk < phislice[etabin - 1].trackidx.size(); itrk++)
0167 clusters[nclust].trackidx.push_back(phislice[etabin - 1].trackidx[itrk]);
0168 }
0169
0170 if (my_pt >= nextbin2_pt && nextbin_pt > 0) {
0171 clusters[nclust].pTtot += nextbin_pt;
0172 clusters[nclust].ntracks += phislice[etabin + 1].ntracks;
0173 clusters[nclust].nxtracks += phislice[etabin + 1].nxtracks;
0174 for (unsigned int itrk = 0; itrk < phislice[etabin + 1].trackidx.size(); itrk++)
0175 clusters[nclust].trackidx.push_back(phislice[etabin + 1].trackidx[itrk]);
0176 phislice[etabin + 1].used = true;
0177 }
0178
0179 nclust++;
0180
0181 }
0182
0183
0184 for (int m = 0; m < nclust - 1; ++m) {
0185 if (((clusters[m + 1].eta - clusters[m].eta) < (3 * etaStep_) / 2) &&
0186 (-(3 * etaStep_) / 2 < (clusters[m + 1].eta - clusters[m].eta))) {
0187 if (clusters[m + 1].pTtot > clusters[m].pTtot) {
0188 clusters[m].eta = clusters[m + 1].eta;
0189 }
0190 clusters[m].pTtot += clusters[m + 1].pTtot;
0191 clusters[m].ntracks += clusters[m + 1].ntracks;
0192 clusters[m].nxtracks += clusters[m + 1].nxtracks;
0193 for (unsigned int itrk = 0; itrk < clusters[m + 1].trackidx.size(); itrk++)
0194 clusters[m].trackidx.push_back(clusters[m + 1].trackidx[itrk]);
0195
0196
0197 for (int m1 = m + 1; m1 < nclust - 1; ++m1)
0198 clusters[m1] = clusters[m1 + 1];
0199
0200 clusters.erase(clusters.begin() + nclust);
0201 nclust--;
0202 }
0203 }
0204
0205 return clusters;
0206 }
0207
0208
0209 template <typename T, typename Pt>
0210 inline void Fill_L2Cluster(T &bin, Pt pt, int ntrk, int ndtrk, std::vector<unsigned int> trkidx) {
0211 bin.pTtot += pt;
0212 bin.ntracks += ntrk;
0213 bin.nxtracks += ndtrk;
0214 for (unsigned int itrk = 0; itrk < trkidx.size(); itrk++)
0215 bin.trackidx.push_back(trkidx[itrk]);
0216 }
0217
0218 inline glbphi_intern DPhi(glbphi_intern phi1, glbphi_intern phi2) {
0219 glbphi_intern x = phi1 - phi2;
0220 if (x >= DoubleToBit(
0221 M_PI, TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, TTTrack_TrackWord::stepPhi0))
0222 x -= DoubleToBit(
0223 2 * M_PI, TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, TTTrack_TrackWord::stepPhi0);
0224 if (x < DoubleToBit(-1 * M_PI,
0225 TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit,
0226 TTTrack_TrackWord::stepPhi0))
0227 x += DoubleToBit(
0228 2 * M_PI, TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, TTTrack_TrackWord::stepPhi0);
0229 return x;
0230 }
0231
0232 inline float DPhi(float phi1, float phi2) {
0233 float x = phi1 - phi2;
0234 if (x >= M_PI)
0235 x -= 2 * M_PI;
0236 if (x < -1 * M_PI)
0237 x += 2 * M_PI;
0238 return x;
0239 }
0240
0241
0242 template <typename T, typename Pt, typename Eta, typename Phi>
0243 inline std::vector<T> L2_clustering(std::vector<std::vector<T> > &L1clusters,
0244 int phiBins_,
0245 Phi phiStep_,
0246 Eta etaStep_) {
0247 std::vector<T> clusters;
0248 for (int phibin = 0; phibin < phiBins_; ++phibin) {
0249 if (L1clusters[phibin].empty())
0250 continue;
0251
0252
0253 sort(L1clusters[phibin].begin(), L1clusters[phibin].end(), [](T &a, T &b) { return a.pTtot > b.pTtot; });
0254
0255 for (unsigned int imax = 0; imax < L1clusters[phibin].size(); ++imax) {
0256 if (L1clusters[phibin][imax].used)
0257 continue;
0258 Pt pt_current = L1clusters[phibin][imax].pTtot;
0259 Pt pt_next = 0;
0260 Pt pt_next2 = 0;
0261 int trk1 = 0;
0262 int trk2 = 0;
0263 int tdtrk1 = 0;
0264 int tdtrk2 = 0;
0265 std::vector<unsigned int> trkidx1;
0266 std::vector<unsigned int> trkidx2;
0267 clusters.push_back(L1clusters[phibin][imax]);
0268
0269 L1clusters[phibin][imax].used = true;
0270
0271
0272 if (phibin == phiBins_ - 1)
0273 continue;
0274
0275 std::vector<unsigned int> used_already;
0276 for (unsigned int icluster = 0; icluster < L1clusters[phibin + 1].size(); ++icluster) {
0277 if (L1clusters[phibin + 1][icluster].used)
0278 continue;
0279
0280 if (((L1clusters[phibin + 1][icluster].eta - L1clusters[phibin][imax].eta) > (3 * etaStep_) / 2) ||
0281 ((L1clusters[phibin + 1][icluster].eta - L1clusters[phibin][imax].eta) < -(3 * etaStep_) / 2))
0282 continue;
0283
0284 pt_next += L1clusters[phibin + 1][icluster].pTtot;
0285 trk1 += L1clusters[phibin + 1][icluster].ntracks;
0286 tdtrk1 += L1clusters[phibin + 1][icluster].nxtracks;
0287
0288 for (unsigned int itrk = 0; itrk < L1clusters[phibin + 1][icluster].trackidx.size(); itrk++)
0289 trkidx1.push_back(L1clusters[phibin + 1][icluster].trackidx[itrk]);
0290 used_already.push_back(icluster);
0291 }
0292
0293 if (pt_next < pt_current) {
0294 Fill_L2Cluster<T, Pt>(clusters[clusters.size() - 1], pt_next, trk1, tdtrk1, trkidx1);
0295 for (unsigned int iused : used_already)
0296 L1clusters[phibin + 1][iused].used = true;
0297 continue;
0298 }
0299
0300 if (phibin == phiBins_ - 2) {
0301 Fill_L2Cluster<T, Pt>(clusters[clusters.size() - 1], pt_next, trk1, tdtrk1, trkidx1);
0302 clusters[clusters.size() - 1].phi = L1clusters[phibin + 1][used_already[0]].phi;
0303 for (unsigned int iused : used_already)
0304 L1clusters[phibin + 1][iused].used = true;
0305 continue;
0306 }
0307 std::vector<int> used_already2;
0308 for (unsigned int icluster = 0; icluster < L1clusters[phibin + 2].size(); ++icluster) {
0309 if (L1clusters[phibin + 2][icluster].used)
0310 continue;
0311 if (((L1clusters[phibin + 2][icluster].eta - L1clusters[phibin][imax].eta) > (3 * etaStep_) / 2) ||
0312 ((L1clusters[phibin + 2][icluster].eta - L1clusters[phibin][imax].eta) < -(3 * etaStep_) / 2))
0313 continue;
0314 pt_next2 += L1clusters[phibin + 2][icluster].pTtot;
0315 trk2 += L1clusters[phibin + 2][icluster].ntracks;
0316 tdtrk2 += L1clusters[phibin + 2][icluster].nxtracks;
0317
0318 for (unsigned int itrk = 0; itrk < L1clusters[phibin + 2][icluster].trackidx.size(); itrk++)
0319 trkidx2.push_back(L1clusters[phibin + 2][icluster].trackidx[itrk]);
0320 used_already2.push_back(icluster);
0321 }
0322 if (pt_next2 < pt_next) {
0323 std::vector<unsigned int> trkidx_both;
0324 trkidx_both.reserve(trkidx1.size() + trkidx2.size());
0325 trkidx_both.insert(trkidx_both.end(), trkidx1.begin(), trkidx1.end());
0326 trkidx_both.insert(trkidx_both.end(), trkidx2.begin(), trkidx2.end());
0327 Fill_L2Cluster<T, Pt>(
0328 clusters[clusters.size() - 1], pt_next + pt_next2, trk1 + trk2, tdtrk1 + tdtrk2, trkidx_both);
0329 clusters[clusters.size() - 1].phi = L1clusters[phibin + 1][used_already[0]].phi;
0330 for (unsigned int iused : used_already)
0331 L1clusters[phibin + 1][iused].used = true;
0332 for (unsigned int iused : used_already2)
0333 L1clusters[phibin + 2][iused].used = true;
0334 }
0335 }
0336 }
0337
0338 int nclust = clusters.size();
0339
0340
0341 for (int m = 0; m < nclust - 1; ++m) {
0342 for (int n = m + 1; n < nclust; ++n) {
0343 if (clusters[n].eta != clusters[m].eta)
0344 continue;
0345 if ((DPhi(clusters[n].phi, clusters[m].phi) > (3 * phiStep_) / 2) ||
0346 (DPhi(clusters[n].phi, clusters[m].phi) < -(3 * phiStep_) / 2))
0347 continue;
0348 if (clusters[n].pTtot > clusters[m].pTtot)
0349 clusters[m].phi = clusters[n].phi;
0350 clusters[m].pTtot += clusters[n].pTtot;
0351 clusters[m].ntracks += clusters[n].ntracks;
0352 clusters[m].nxtracks += clusters[n].nxtracks;
0353 for (unsigned int itrk = 0; itrk < clusters[n].trackidx.size(); itrk++)
0354 clusters[m].trackidx.push_back(clusters[n].trackidx[itrk]);
0355 for (int m1 = n; m1 < nclust - 1; ++m1)
0356 clusters[m1] = clusters[m1 + 1];
0357 clusters.erase(clusters.begin() + nclust);
0358
0359 nclust--;
0360 }
0361 }
0362 return clusters;
0363 }
0364 }
0365 #endif