L1TTrackMatch/interface/L1Clustering.h

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 //Each individual box in the eta and phi dimension.
0009 //  Also used to store final cluster data for each zbin.
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;  //average phi value (halfway b/t min and max)
0018   float eta;  //average eta value
0019   std::vector<unsigned int> trackidx;
0020 };
0021
0022 //store important information for plots
0023 struct MaxZBin {
0024   int znum;    //Numbered from 0 to nzbins (16, 32, or 64) in order
0025   int nclust;  //number of clusters in this bin
0026   float zbincenter;
0027   std::vector<EtaPhiBin> clusters;  //list of all the clusters in this bin
0028   float ht;                         //sum of all cluster pTs--only the zbin with the maximum ht is stored
0029 };
0030
0031 inline std::vector<EtaPhiBin> L1_clustering(EtaPhiBin *phislice, int etaBins_, float etaStep_) {
0032   std::vector<EtaPhiBin> clusters;
0033   // Find eta-phibin with maxpT, make center of cluster, add neighbors if not already used
0034   int nclust = 0;
0035
0036   // get tracks in eta bins in increasing eta order
0037   for (int etabin = 0; etabin < etaBins_; ++etabin) {
0038     float my_pt = 0, previousbin_pt = 0;  //, nextbin_pt=0, next2bin_pt=0;
0039     float nextbin_pt = 0, nextbin2_pt = 0;
0040
0041     // skip (already) used tracks
0042     if (phislice[etabin].used)
0043       continue;
0044     my_pt = phislice[etabin].pTtot;
0045     if (my_pt == 0)
0046       continue;
0047     //get previous bin pT
0048     if (etabin > 0 && !phislice[etabin - 1].used)
0049       previousbin_pt = phislice[etabin - 1].pTtot;
0050
0051     // get next bins pt
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     // check if pT of current cluster is higher than neighbors
0059     if (my_pt < previousbin_pt || my_pt <= nextbin_pt) {
0060       // if unused pT in the left neighbor, spit it out as a cluster
0061       if (previousbin_pt > 0) {
0062         clusters.push_back(phislice[etabin - 1]);
0063         phislice[etabin - 1].used = true;
0064         nclust++;
0065       }
0066       continue;  //if it is not the local max pT skip
0067     }
0068     // here reach only unused local max clusters
0069     clusters.push_back(phislice[etabin]);
0070     phislice[etabin].used = true;  //if current bin a cluster
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   }  // for each etabin
0091
0092   // Merge close-by clusters
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;  // total ntrk
0100       clusters[m].numtdtrks += clusters[m + 1].numtdtrks;  // total ndisp
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       // if remove the merged cluster - all the others must be closer to 0
0105       for (int m1 = m + 1; m1 < nclust - 1; ++m1) {
0106         clusters[m1] = clusters[m1 + 1];
0107         //clusters.erase(clusters.begin()+m1);
0108       }
0109       //  clusters[m1] = clusters[m1 + 1];
0110       clusters.erase(clusters.begin() + nclust);
0111       nclust--;
0112       m = -1;
0113     }  // end if clusters neighbor in eta
0114   }    // end for (m) loop
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) {  //Find eta-phibin with highest pT
0142     if (L1clusters[phibin].empty())
0143       continue;
0144
0145     // sort L1 clusters max -> min
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;  //current cluster (pt0)
0153       float pt_next = 0;                                  // next phi bin (pt1)
0154       float pt_next2 = 0;                                 // next to next phi bin2 (pt2)
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       // if we are in the last phi bin, dont check phi+1 phi+2
0166       if (phibin == phiBins_ - 1)
0167         continue;
0168       std::vector<unsigned int> used_already;  //keep phi+1 clusters that have been used
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) {  // if pt1<pt1, merge both clusters
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       // if phi = next to last bin there is no "next to next"
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;  //keep used clusters in phi+2
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     }  // for each L1 cluster
0222   }    // for each phibin
0223
0224   int nclust = clusters.size();
0225
0226   // merge close-by clusters
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     }  // end of n-loop
0248   }    // end of m-loop
0249   return clusters;
0250 }
0251 #endif