File indexing completed on 2022-12-25 23:44:18
0001 #ifndef L1Trigger_L1TTrackMatch_L1Clustering_HH
0002 #define L1Trigger_L1TTrackMatch_L1Clustering_HH
0003 #include <cmath>
0004 #include <cstdlib>
0005 #include <vector>
0006 #include <algorithm>
0007
0008
0009
0010 struct EtaPhiBin {
0011 float pTtot = 0;
0012 int numtracks = 0;
0013 int numttrks = 0;
0014 int numtdtrks = 0;
0015 int numttdtrks = 0;
0016 bool used = false;
0017 float phi;
0018 float eta;
0019 std::vector<unsigned int> trackidx;
0020 };
0021
0022
0023 struct MaxZBin {
0024 int znum;
0025 int nclust;
0026 float zbincenter;
0027 std::vector<EtaPhiBin> clusters;
0028 float ht;
0029 };
0030
0031 inline std::vector<EtaPhiBin> L1_clustering(EtaPhiBin *phislice, int etaBins_, float etaStep_) {
0032 std::vector<EtaPhiBin> clusters;
0033
0034 int nclust = 0;
0035
0036
0037 for (int etabin = 0; etabin < etaBins_; ++etabin) {
0038 float my_pt = 0, previousbin_pt = 0;
0039 float nextbin_pt = 0, nextbin2_pt = 0;
0040
0041
0042 if (phislice[etabin].used)
0043 continue;
0044 my_pt = phislice[etabin].pTtot;
0045 if (my_pt == 0)
0046 continue;
0047
0048 if (etabin > 0 && !phislice[etabin - 1].used)
0049 previousbin_pt = phislice[etabin - 1].pTtot;
0050
0051
0052 if (etabin < etaBins_ - 1 && !phislice[etabin + 1].used) {
0053 nextbin_pt = phislice[etabin + 1].pTtot;
0054 if (etabin < etaBins_ - 2 && !phislice[etabin + 2].used) {
0055 nextbin2_pt = phislice[etabin + 2].pTtot;
0056 }
0057 }
0058
0059 if (my_pt < previousbin_pt || my_pt <= nextbin_pt) {
0060
0061 if (previousbin_pt > 0) {
0062 clusters.push_back(phislice[etabin - 1]);
0063 phislice[etabin - 1].used = true;
0064 nclust++;
0065 }
0066 continue;
0067 }
0068
0069 clusters.push_back(phislice[etabin]);
0070 phislice[etabin].used = true;
0071 if (previousbin_pt > 0) {
0072 clusters[nclust].pTtot += previousbin_pt;
0073 clusters[nclust].numtracks += phislice[etabin - 1].numtracks;
0074 clusters[nclust].numtdtrks += phislice[etabin - 1].numtdtrks;
0075 for (unsigned int itrk = 0; itrk < phislice[etabin - 1].trackidx.size(); itrk++)
0076 clusters[nclust].trackidx.push_back(phislice[etabin - 1].trackidx[itrk]);
0077 }
0078
0079 if (my_pt >= nextbin2_pt && nextbin_pt > 0) {
0080 clusters[nclust].pTtot += nextbin_pt;
0081 clusters[nclust].numtracks += phislice[etabin + 1].numtracks;
0082 clusters[nclust].numtdtrks += phislice[etabin + 1].numtdtrks;
0083 for (unsigned int itrk = 0; itrk < phislice[etabin + 1].trackidx.size(); itrk++)
0084 clusters[nclust].trackidx.push_back(phislice[etabin + 1].trackidx[itrk]);
0085 phislice[etabin + 1].used = true;
0086 }
0087
0088 nclust++;
0089
0090 }
0091
0092
0093 for (int m = 0; m < nclust - 1; ++m) {
0094 if (std::abs(clusters[m + 1].eta - clusters[m].eta) < 1.5 * etaStep_) {
0095 if (clusters[m + 1].pTtot > clusters[m].pTtot) {
0096 clusters[m].eta = clusters[m + 1].eta;
0097 }
0098 clusters[m].pTtot += clusters[m + 1].pTtot;
0099 clusters[m].numtracks += clusters[m + 1].numtracks;
0100 clusters[m].numtdtrks += clusters[m + 1].numtdtrks;
0101 for (unsigned int itrk = 0; itrk < clusters[m + 1].trackidx.size(); itrk++)
0102 clusters[m].trackidx.push_back(clusters[m + 1].trackidx[itrk]);
0103
0104
0105 for (int m1 = m + 1; m1 < nclust - 1; ++m1) {
0106 clusters[m1] = clusters[m1 + 1];
0107
0108 }
0109
0110 clusters.erase(clusters.begin() + nclust);
0111 nclust--;
0112 m = -1;
0113 }
0114 }
0115
0116 return clusters;
0117 }
0118
0119 inline void Fill_L2Cluster(EtaPhiBin &bin, float pt, int ntrk, int ndtrk, std::vector<unsigned int> trkidx) {
0120 bin.pTtot += pt;
0121 bin.numtracks += ntrk;
0122 bin.numtdtrks += ndtrk;
0123 for (unsigned int itrk = 0; itrk < trkidx.size(); itrk++)
0124 bin.trackidx.push_back(trkidx[itrk]);
0125 }
0126
0127 inline float DPhi(float phi1, float phi2) {
0128 float x = phi1 - phi2;
0129 if (x >= M_PI)
0130 x -= 2 * M_PI;
0131 if (x < -1 * M_PI)
0132 x += 2 * M_PI;
0133 return x;
0134 }
0135
0136 inline std::vector<EtaPhiBin> L2_clustering(std::vector<std::vector<EtaPhiBin>> &L1clusters,
0137 int phiBins_,
0138 float phiStep_,
0139 float etaStep_) {
0140 std::vector<EtaPhiBin> clusters;
0141 for (int phibin = 0; phibin < phiBins_; ++phibin) {
0142 if (L1clusters[phibin].empty())
0143 continue;
0144
0145
0146 sort(L1clusters[phibin].begin(), L1clusters[phibin].end(), [](struct EtaPhiBin &a, struct EtaPhiBin &b) {
0147 return a.pTtot > b.pTtot;
0148 });
0149 for (unsigned int imax = 0; imax < L1clusters[phibin].size(); ++imax) {
0150 if (L1clusters[phibin][imax].used)
0151 continue;
0152 float pt_current = L1clusters[phibin][imax].pTtot;
0153 float pt_next = 0;
0154 float pt_next2 = 0;
0155 int trk1 = 0;
0156 int trk2 = 0;
0157 int tdtrk1 = 0;
0158 int tdtrk2 = 0;
0159 std::vector<unsigned int> trkidx1;
0160 std::vector<unsigned int> trkidx2;
0161 clusters.push_back(L1clusters[phibin][imax]);
0162
0163 L1clusters[phibin][imax].used = true;
0164
0165
0166 if (phibin == phiBins_ - 1)
0167 continue;
0168 std::vector<unsigned int> used_already;
0169 for (unsigned int icluster = 0; icluster < L1clusters[phibin + 1].size(); ++icluster) {
0170 if (L1clusters[phibin + 1][icluster].used)
0171 continue;
0172 if (std::abs(L1clusters[phibin + 1][icluster].eta - L1clusters[phibin][imax].eta) > 1.5 * etaStep_)
0173 continue;
0174 pt_next += L1clusters[phibin + 1][icluster].pTtot;
0175 trk1 += L1clusters[phibin + 1][icluster].numtracks;
0176 tdtrk1 += L1clusters[phibin + 1][icluster].numtdtrks;
0177 for (unsigned int itrk = 0; itrk < L1clusters[phibin + 1][icluster].trackidx.size(); itrk++)
0178 trkidx1.push_back(L1clusters[phibin + 1][icluster].trackidx[itrk]);
0179 used_already.push_back(icluster);
0180 }
0181
0182 if (pt_next < pt_current) {
0183 Fill_L2Cluster(clusters[clusters.size() - 1], pt_next, trk1, tdtrk1, trkidx1);
0184 for (unsigned int iused : used_already)
0185 L1clusters[phibin + 1][iused].used = true;
0186 continue;
0187 }
0188
0189 if (phibin == phiBins_ - 2) {
0190 Fill_L2Cluster(clusters[clusters.size() - 1], pt_next, trk1, tdtrk1, trkidx1);
0191 clusters[clusters.size() - 1].phi = L1clusters[phibin + 1][used_already[0]].phi;
0192 for (unsigned int iused : used_already)
0193 L1clusters[phibin + 1][iused].used = true;
0194 continue;
0195 }
0196 std::vector<int> used_already2;
0197 for (unsigned int icluster = 0; icluster < L1clusters[phibin + 2].size(); ++icluster) {
0198 if (L1clusters[phibin + 2][icluster].used)
0199 continue;
0200 if (std::abs(L1clusters[phibin + 2][icluster].eta - L1clusters[phibin][imax].eta) > 1.5 * etaStep_)
0201 continue;
0202 pt_next2 += L1clusters[phibin + 2][icluster].pTtot;
0203 trk2 += L1clusters[phibin + 2][icluster].numtracks;
0204 tdtrk2 += L1clusters[phibin + 2][icluster].numtdtrks;
0205 for (unsigned int itrk = 0; itrk < L1clusters[phibin + 2][icluster].trackidx.size(); itrk++)
0206 trkidx2.push_back(L1clusters[phibin + 2][icluster].trackidx[itrk]);
0207 used_already2.push_back(icluster);
0208 }
0209 if (pt_next2 < pt_next) {
0210 std::vector<unsigned int> trkidx_both;
0211 trkidx_both.reserve(trkidx1.size() + trkidx2.size());
0212 trkidx_both.insert(trkidx_both.end(), trkidx1.begin(), trkidx1.end());
0213 trkidx_both.insert(trkidx_both.end(), trkidx2.begin(), trkidx2.end());
0214 Fill_L2Cluster(clusters[clusters.size() - 1], pt_next + pt_next2, trk1 + trk2, tdtrk1 + tdtrk2, trkidx_both);
0215 clusters[clusters.size() - 1].phi = L1clusters[phibin + 1][used_already[0]].phi;
0216 for (unsigned int iused : used_already)
0217 L1clusters[phibin + 1][iused].used = true;
0218 for (unsigned int iused : used_already2)
0219 L1clusters[phibin + 2][iused].used = true;
0220 }
0221 }
0222 }
0223
0224 int nclust = clusters.size();
0225
0226
0227 for (int m = 0; m < nclust - 1; ++m) {
0228 for (int n = m + 1; n < nclust; ++n) {
0229 if (clusters[n].eta != clusters[m].eta)
0230 continue;
0231 if (std::abs(DPhi(clusters[n].phi, clusters[m].phi)) > 1.5 * phiStep_)
0232 continue;
0233
0234 if (clusters[n].pTtot > clusters[m].pTtot)
0235 clusters[m].phi = clusters[n].phi;
0236
0237 clusters[m].pTtot += clusters[n].pTtot;
0238 clusters[m].numtracks += clusters[n].numtracks;
0239 clusters[m].numtdtrks += clusters[n].numtdtrks;
0240 for (unsigned int itrk = 0; itrk < clusters[n].trackidx.size(); itrk++)
0241 clusters[m].trackidx.push_back(clusters[n].trackidx[itrk]);
0242 for (int m1 = n; m1 < nclust - 1; ++m1)
0243 clusters[m1] = clusters[m1 + 1];
0244 clusters.erase(clusters.begin() + nclust);
0245
0246 nclust--;
0247 }
0248 }
0249 return clusters;
0250 }
0251 #endif