L1TTrackMatch/plugins/L1TrackJetClustering.h

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   //For precision studies
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;  // 40cm / 0.1
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;  // The remove 1 bit from nBits to account for the sign
0029     if (digitized_value > digitized_maximum)
0030       digitized_value = digitized_maximum;
0031     if (value < 0)
0032       digitized_value = (1 << maxBits) - digitized_value;  // two's complement encoding
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))) {  // check the sign
0039       isign = -1;
0040       bits = (1 << (maxBits + 1)) - bits;
0041     }
0042     return (double(bits & digitized_maximum) + 0.5) * step * isign;
0043   }
0044
0045   // eta/phi clusters - simulation
0046   struct EtaPhiBin {
0047     float pTtot;
0048     int ntracks;
0049     int nxtracks;
0050     bool used;
0051     float phi;  //average phi value (halfway b/t min and max)
0052     float eta;  //average eta value
0053     std::vector<unsigned int> trackidx;
0054   };
0055   // z bin struct - simulation (used if z bin are many)
0056   struct MaxZBin {
0057     int znum;    //Numbered from 0 to nzbins (16, 32, or 64) in order
0058     int nclust;  //number of clusters in this bin
0059     float zbincenter;
0060     std::vector<EtaPhiBin> clusters;  //list of all the clusters in this bin
0061     float ht;                         //sum of all cluster pTs--only the zbin with the maximum ht is stored
0062   };
0063
0064   // eta/phi clusters - emulation
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;  //average phi value (halfway b/t min and max)
0071     glbeta_intern eta;  //average eta value
0072     std::vector<unsigned int> trackidx;
0073   };
0074
0075   // z bin struct - emulation (used if z bin are many)
0076   struct TrackJetEmulationMaxZBin {
0077     int znum;    //Numbered from 0 to nzbins (16, 32, or 64) in order
0078     int nclust;  //number of clusters in this bin
0079     z0_intern zbincenter;
0080     std::vector<TrackJetEmulationEtaPhiBin> clusters;  //list of all the clusters in this bin
0081     pt_intern ht;  //sum of all cluster pTs--only the zbin with the maximum ht is stored
0082   };
0083
0084   // track quality cuts
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_)  // 4 stubs are the lowest track quality and have different cuts
0101         PassQuality = true;
0102       if (trk_nstub > 4 && trk_bendchi2 < nStubs5PromptBend_ &&
0103           trk_chi2 < nStubs5PromptChi2_)  // above 4 stubs diffent selection imposed (genrally looser)
0104         PassQuality = true;
0105     } else {
0106       if (trk_nstub == 4 && trk_bendchi2 < nStubs4DisplacedBend_ &&
0107           trk_chi2 < nStubs4DisplacedChi2_)  // 4 stubs are the lowest track quality and have different cuts
0108         PassQuality = true;
0109       if (trk_nstub > 4 && trk_bendchi2 < nStubs5DisplacedBend_ &&
0110           trk_chi2 < nStubs5DisplacedChi2_)  // above 4 stubs diffent selection imposed (genrally looser)
0111         PassQuality = true;
0112     }
0113     return PassQuality;
0114   }
0115
0116   // L1 clustering (in eta)
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     // Find eta bin with maxpT, make center of cluster, add neighbors if not already used
0121     int nclust = 0;
0122
0123     // get tracks in eta bins in increasing eta order
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       // skip (already) used tracks
0131       if (phislice[etabin].used)
0132         continue;
0133
0134       my_pt = phislice[etabin].pTtot;
0135       if (my_pt == 0)
0136         continue;
0137
0138       //get previous bin pT
0139       if (etabin > 0 && !phislice[etabin - 1].used)
0140         previousbin_pt = phislice[etabin - 1].pTtot;
0141
0142       // get next bins pt
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       // check if pT of current cluster is higher than neighbors
0150       if (my_pt < previousbin_pt || my_pt <= nextbin_pt) {
0151         // if unused pT in the left neighbor, spit it out as a cluster
0152         if (previousbin_pt > 0) {
0153           clusters.push_back(phislice[etabin - 1]);
0154           phislice[etabin - 1].used = true;
0155           nclust++;
0156         }
0157         continue;  //if it is not the local max pT skip
0158       }
0159       // here reach only unused local max clusters
0160       clusters.push_back(phislice[etabin]);
0161       phislice[etabin].used = true;  //if current bin a cluster
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     }  // for each etabin
0182
0183     // Merge close-by clusters
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;    // total ntrk
0192         clusters[m].nxtracks += clusters[m + 1].nxtracks;  // total ndisp
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         // if remove the merged cluster - all the others must be closer to 0
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       }  // end if for cluster merging
0203     }    // end for (m) loop
0204
0205     return clusters;
0206   }
0207
0208   // Fill L2 clusters (helper function)
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   // L2 clustering (in phi)
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) {  //Find eta-phibin with highest pT
0249       if (L1clusters[phibin].empty())
0250         continue;
0251
0252       // sort L1 clusters max -> min
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;  //current cluster (pt0)
0259         Pt pt_next = 0;                                  // next phi bin (pt1)
0260         Pt pt_next2 = 0;                                 // next to next phi bin2 (pt2)
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         // if we are in the last phi bin, dont check phi+1 phi+2
0272         if (phibin == phiBins_ - 1)
0273           continue;
0274
0275         std::vector<unsigned int> used_already;  //keep phi+1 clusters that have been used
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) {  // if pt1<pt1, merge both clusters
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         // if phi = next to last bin there is no "next to next"
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;  //keep used clusters in phi+2
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       }  // for each L1 cluster
0336     }    // for each phibin
0337
0338     int nclust = clusters.size();
0339
0340     // merge close-by clusters
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       }  // end of n-loop
0361     }    // end of m-loop
0362     return clusters;
0363   }
0364 }  // namespace l1ttrackjet
0365 #endif