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   static constexpr int kEtaWordLength = 15;
0021   static constexpr int kPhiWordLength = 12;
0022
0023   //Constants used for jet clustering in eta direction
0024   static constexpr int kThirteenBitMask = 0b1111111111111;
0025   static constexpr int kEtaFineBinEdge1 = 0b0011001100110;
0026   static constexpr int kEtaFineBinEdge2 = 0b0110011001100;
0027   static constexpr int kEtaFineBinEdge3 = 0b1001100110010;
0028   static constexpr int kEtaFineBinEdge4 = 0b1100110011000;
0029
0030   //Constants used for jet clustering in phi direction
0031   static constexpr int kTwelveBitMask = 0b011111111111;
0032   static constexpr int kPhiBinHalfWidth = 0b000100101111;
0033   static constexpr int kNumPhiBins = 27;
0034   static constexpr int kPhiBinZeroOffset = 12;  // phi bin zero offset between firmware and emulator
0035
0036   typedef ap_ufixed<TTTrack_TrackWord::TrackBitWidths::kRinvSize - 1, PT_INTPART_BITS, AP_TRN, AP_SAT> pt_intern;
0037   typedef ap_fixed<TTTrack_TrackWord::TrackBitWidths::kTanlSize, ETA_INTPART_BITS, AP_TRN, AP_SAT> glbeta_intern;
0038   typedef ap_int<TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit> glbphi_intern;
0039   typedef ap_int<TTTrack_TrackWord::TrackBitWidths::kZ0Size> z0_intern;  // 40cm / 0.1
0040   typedef ap_uint<TTTrack_TrackWord::TrackBitWidths::kD0Size> d0_intern;
0041
0042   inline const unsigned int DoubleToBit(double value, unsigned int maxBits, double step) {
0043     unsigned int digitized_value = std::floor(std::abs(value) / step);
0044     unsigned int digitized_maximum = (1 << (maxBits - 1)) - 1;  // The remove 1 bit from nBits to account for the sign
0045     if (digitized_value > digitized_maximum)
0046       digitized_value = digitized_maximum;
0047     if (value < 0)
0048       digitized_value = (1 << maxBits) - digitized_value;  // two's complement encoding
0049     return digitized_value;
0050   }
0051   inline const double BitToDouble(unsigned int bits, unsigned int maxBits, double step) {
0052     int isign = 1;
0053     unsigned int digitized_maximum = (1 << maxBits) - 1;
0054     if (bits & (1 << (maxBits - 1))) {  // check the sign
0055       isign = -1;
0056       bits = (1 << (maxBits + 1)) - bits;
0057     }
0058     return (double(bits & digitized_maximum) + 0.5) * step * isign;
0059   }
0060
0061   // eta/phi clusters - simulation
0062   struct EtaPhiBin {
0063     float pTtot;
0064     int ntracks;
0065     int nxtracks;
0066     bool used;
0067     float phi;  //average phi value (halfway b/t min and max)
0068     float eta;  //average eta value
0069     std::vector<unsigned int> trackidx;
0070   };
0071   // z bin struct - simulation (used if z bin are many)
0072   struct MaxZBin {
0073     int znum;    //Numbered from 0 to nzbins (16, 32, or 64) in order
0074     int nclust;  //number of clusters in this bin
0075     float zbincenter;
0076     std::vector<EtaPhiBin> clusters;  //list of all the clusters in this bin
0077     float ht;                         //sum of all cluster pTs--only the zbin with the maximum ht is stored
0078   };
0079
0080   // eta/phi clusters - emulation
0081   struct TrackJetEmulationEtaPhiBin {
0082     pt_intern pTtot;
0083     l1t::TkJetWord::nt_t ntracks;
0084     l1t::TkJetWord::nx_t nxtracks;
0085     bool used;
0086     glbphi_intern phi;  //average phi value (halfway b/t min and max)
0087     glbeta_intern eta;  //average eta value
0088     std::vector<unsigned int> trackidx;
0089   };
0090
0091   // z bin struct - emulation (used if z bin are many)
0092   struct TrackJetEmulationMaxZBin {
0093     int znum;    //Numbered from 0 to nzbins (16, 32, or 64) in order
0094     int nclust;  //number of clusters in this bin
0095     z0_intern zbincenter;
0096     std::vector<TrackJetEmulationEtaPhiBin> clusters;  //list of all the clusters in this bin
0097     pt_intern ht;  //sum of all cluster pTs--only the zbin with the maximum ht is stored
0098   };
0099
0100   // track quality cuts
0101   inline bool TrackQualitySelection(int trk_nstub,
0102                                     double trk_chi2,
0103                                     double trk_bendchi2,
0104                                     double nStubs4PromptBend_,
0105                                     double nStubs5PromptBend_,
0106                                     double nStubs4PromptChi2_,
0107                                     double nStubs5PromptChi2_,
0108                                     double nStubs4DisplacedBend_,
0109                                     double nStubs5DisplacedBend_,
0110                                     double nStubs4DisplacedChi2_,
0111                                     double nStubs5DisplacedChi2_,
0112                                     bool displaced_) {
0113     bool PassQuality = false;
0114     if (!displaced_) {
0115       if (trk_nstub == 4 && trk_bendchi2 < nStubs4PromptBend_ &&
0116           trk_chi2 < nStubs4PromptChi2_)  // 4 stubs are the lowest track quality and have different cuts
0117         PassQuality = true;
0118       if (trk_nstub > 4 && trk_bendchi2 < nStubs5PromptBend_ &&
0119           trk_chi2 < nStubs5PromptChi2_)  // above 4 stubs diffent selection imposed (genrally looser)
0120         PassQuality = true;
0121     } else {
0122       if (trk_nstub == 4 && trk_bendchi2 < nStubs4DisplacedBend_ &&
0123           trk_chi2 < nStubs4DisplacedChi2_)  // 4 stubs are the lowest track quality and have different cuts
0124         PassQuality = true;
0125       if (trk_nstub > 4 && trk_bendchi2 < nStubs5DisplacedBend_ &&
0126           trk_chi2 < nStubs5DisplacedChi2_)  // above 4 stubs diffent selection imposed (genrally looser)
0127         PassQuality = true;
0128     }
0129     return PassQuality;
0130   }
0131
0132   // Eta binning following the hardware logic
0133   inline unsigned int eta_bin_firmwareStyle(int eta_word) {
0134     //Function that reads in 15-bit eta word (in two's complement) and returns the index of eta bin in which it belongs
0135     //Logic follows exactly according to the firmware
0136     //We will first determine if eta is pos/neg from the first bit. Each half of the grid is then split into four coarse bins. Bits 2&3 determine which coarse bin to assign
0137     //The coarse bins are split into 5 fine bins, final 13 bits determine which of these coarse bins this track needs to assign
0138     int eta_coarse = 0;
0139     int eta_fine = 0;
0140
0141     if (eta_word & (1 << kEtaWordLength)) {  //If eta is negative (first/leftmost bit is 1)
0142       //Second and third bits contain information about which of four coarse bins
0143       eta_coarse = 5 * ((eta_word & (3 << (kEtaWordLength - 2))) >> (kEtaWordLength - 2));
0144     } else {  //else eta is positive (first/leftmost bit is 0)
0145       eta_coarse = 5 * (4 + ((eta_word & (3 << (kEtaWordLength - 2))) >> (kEtaWordLength - 2)));
0146     }
0147
0148     //Now get the fine bin index. The numbers correspond to the decimal representation of fine bin edges in binary
0149     int j = eta_word & kThirteenBitMask;
0150     if (j < kEtaFineBinEdge1)
0151       eta_fine = 0;
0152     else if (j < kEtaFineBinEdge2)
0153       eta_fine = 1;
0154     else if (j < kEtaFineBinEdge3)
0155       eta_fine = 2;
0156     else if (j < kEtaFineBinEdge4)
0157       eta_fine = 3;
0158     else
0159       eta_fine = 4;
0160
0161     //Final eta bin is coarse bin index + fine bin index, subtract 8 to make eta_bin at eta=-2.4 have index=0
0162     int eta_ = (eta_coarse + eta_fine) - 8;
0163     return eta_;
0164   }
0165
0166   // Phi binning following the hardware logic
0167   inline unsigned int phi_bin_firmwareStyle(int phi_sector_raw, int phi_word) {
0168     //Function that reads in decimal integers phi_sector_raw and phi_word and returns the index of phi bin in which it belongs
0169     //phi_sector_raw is integer 0-8 correspoding to one of 9 phi sectors
0170     //phi_word is 12 bit word representing the phi value measured w.r.t the center of the sector
0171
0172     int phi_coarse = 3 * phi_sector_raw;  //phi_coarse is index of phi coarse binning (sector edges)
0173     int phi_fine = 0;  //phi_fine is index of fine bin inside sectors. Each sector contains 3 fine bins
0174
0175     //Determine fine bin. First bit is sign, next 11 bits determine fine bin
0176     //303 is distance from 0 to first fine bin edge
0177     //2047 is eleven 1's, use the &2047 to extract leftmost 11 bits.
0178
0179     //The allowed range for phi goes further than the edges of bin 0 or 2 (bit value 909). There's an apparent risk of phi being > 909, however this will always mean the track is in the next link (i.e. track beyond bin 2 in this link means track is actually in bin 0 of adjacent link)
0180
0181     if (phi_word & (1 << (kPhiWordLength - 1))) {  //if phi is negative (first bit 1)
0182       //Since negative, we 'flip' the phi word, then check if it is in fine bin 0 or 1
0183       if ((kTwelveBitMask - (phi_word & kTwelveBitMask)) > kPhiBinHalfWidth) {
0184         phi_fine = 0;
0185       } else if ((kTwelveBitMask - (phi_word & kTwelveBitMask)) < kPhiBinHalfWidth) {
0186         phi_fine = 1;
0187       }
0188     } else {  //else phi is positive (first bit 0)
0189       //positive phi, no 'flip' necessary. Just check if in  fine bin 1 or 2
0190       if ((phi_word & kTwelveBitMask) < kPhiBinHalfWidth) {
0191         phi_fine = 1;
0192       } else if ((phi_word & kTwelveBitMask) > kPhiBinHalfWidth) {
0193         phi_fine = 2;
0194       }
0195     }
0196
0197     // Final operation is a shift by pi (half a grid) to make bin at index=0 at -pi
0198     int phi_bin_ = (phi_coarse + phi_fine + kPhiBinZeroOffset) % kNumPhiBins;
0199
0200     return phi_bin_;
0201   }
0202
0203   // L1 clustering (in eta)
0204   template <typename T, typename Pt, typename Eta, typename Phi>
0205   inline std::vector<T> L1_clustering(T *phislice, int etaBins_, Eta etaStep_) {
0206     std::vector<T> clusters;
0207     // Find eta bin with maxpT, make center of cluster, add neighbors if not already used
0208     int nclust = 0;
0209
0210     // get tracks in eta bins in increasing eta order
0211     for (int etabin = 0; etabin < etaBins_; ++etabin) {
0212       Pt my_pt = 0;
0213       Pt previousbin_pt = 0;
0214       Pt nextbin_pt = 0;
0215       Pt nextbin2_pt = 0;
0216
0217       // skip (already) used tracks
0218       if (phislice[etabin].used)
0219         continue;
0220
0221       my_pt = phislice[etabin].pTtot;
0222       if (my_pt == 0)
0223         continue;
0224
0225       //get previous bin pT
0226       if (etabin > 0 && !phislice[etabin - 1].used)
0227         previousbin_pt = phislice[etabin - 1].pTtot;
0228
0229       // get next bins pt
0230       if (etabin < etaBins_ - 1 && !phislice[etabin + 1].used) {
0231         nextbin_pt = phislice[etabin + 1].pTtot;
0232         if (etabin < etaBins_ - 2 && !phislice[etabin + 2].used) {
0233           nextbin2_pt = phislice[etabin + 2].pTtot;
0234         }
0235       }
0236       // check if pT of current cluster is higher than neighbors
0237       if (my_pt < previousbin_pt || my_pt <= nextbin_pt) {
0238         // if unused pT in the left neighbor, spit it out as a cluster
0239         if (previousbin_pt > 0) {
0240           clusters.push_back(phislice[etabin - 1]);
0241           phislice[etabin - 1].used = true;
0242           nclust++;
0243         }
0244         continue;  //if it is not the local max pT skip
0245       }
0246       // here reach only unused local max clusters
0247       clusters.push_back(phislice[etabin]);
0248       phislice[etabin].used = true;  //if current bin a cluster
0249       if (previousbin_pt > 0) {
0250         clusters[nclust].pTtot += previousbin_pt;
0251         clusters[nclust].ntracks += phislice[etabin - 1].ntracks;
0252         clusters[nclust].nxtracks += phislice[etabin - 1].nxtracks;
0253         for (unsigned int itrk = 0; itrk < phislice[etabin - 1].trackidx.size(); itrk++)
0254           clusters[nclust].trackidx.push_back(phislice[etabin - 1].trackidx[itrk]);
0255       }
0256
0257       if (my_pt >= nextbin2_pt && nextbin_pt > 0) {
0258         clusters[nclust].pTtot += nextbin_pt;
0259         clusters[nclust].ntracks += phislice[etabin + 1].ntracks;
0260         clusters[nclust].nxtracks += phislice[etabin + 1].nxtracks;
0261         for (unsigned int itrk = 0; itrk < phislice[etabin + 1].trackidx.size(); itrk++)
0262           clusters[nclust].trackidx.push_back(phislice[etabin + 1].trackidx[itrk]);
0263         phislice[etabin + 1].used = true;
0264       }
0265
0266       nclust++;
0267
0268     }  // for each etabin
0269
0270     // Merge close-by clusters
0271     for (int m = 0; m < nclust - 1; ++m) {
0272       if (((clusters[m + 1].eta - clusters[m].eta) < (3 * etaStep_) / 2) &&
0273           (-(3 * etaStep_) / 2 < (clusters[m + 1].eta - clusters[m].eta))) {
0274         if (clusters[m + 1].pTtot > clusters[m].pTtot) {
0275           clusters[m].eta = clusters[m + 1].eta;
0276         }
0277         clusters[m].pTtot += clusters[m + 1].pTtot;
0278         clusters[m].ntracks += clusters[m + 1].ntracks;    // total ntrk
0279         clusters[m].nxtracks += clusters[m + 1].nxtracks;  // total ndisp
0280         for (unsigned int itrk = 0; itrk < clusters[m + 1].trackidx.size(); itrk++)
0281           clusters[m].trackidx.push_back(clusters[m + 1].trackidx[itrk]);
0282
0283         // if remove the merged cluster - all the others must be closer to 0
0284         for (int m1 = m + 1; m1 < nclust - 1; ++m1)
0285           clusters[m1] = clusters[m1 + 1];
0286
0287         clusters.erase(clusters.begin() + nclust);
0288         nclust--;
0289       }  // end if for cluster merging
0290     }  // end for (m) loop
0291
0292     return clusters;
0293   }
0294
0295   // Fill L2 clusters (helper function)
0296   template <typename T, typename Pt>
0297   inline void Fill_L2Cluster(T &bin, Pt pt, int ntrk, int ndtrk, std::vector<unsigned int> trkidx) {
0298     bin.pTtot += pt;
0299     bin.ntracks += ntrk;
0300     bin.nxtracks += ndtrk;
0301     for (unsigned int itrk = 0; itrk < trkidx.size(); itrk++)
0302       bin.trackidx.push_back(trkidx[itrk]);
0303   }
0304
0305   inline glbphi_intern DPhi(glbphi_intern phi1, glbphi_intern phi2) {
0306     glbphi_intern x = phi1 - phi2;
0307     if (x >= DoubleToBit(
0308                  M_PI, TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, TTTrack_TrackWord::stepPhi0))
0309       x -= DoubleToBit(
0310           2 * M_PI, TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, TTTrack_TrackWord::stepPhi0);
0311     if (x < DoubleToBit(-1 * M_PI,
0312                         TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit,
0313                         TTTrack_TrackWord::stepPhi0))
0314       x += DoubleToBit(
0315           2 * M_PI, TTTrack_TrackWord::TrackBitWidths::kPhiSize + kExtraGlobalPhiBit, TTTrack_TrackWord::stepPhi0);
0316     return x;
0317   }
0318
0319   inline float DPhi(float phi1, float phi2) {
0320     float x = phi1 - phi2;
0321     if (x >= M_PI)
0322       x -= 2 * M_PI;
0323     if (x < -1 * M_PI)
0324       x += 2 * M_PI;
0325     return x;
0326   }
0327
0328   // L2 clustering (in phi)
0329   template <typename T, typename Pt, typename Eta, typename Phi>
0330   inline std::vector<T> L2_clustering(std::vector<std::vector<T> > &L1clusters,
0331                                       int phiBins_,
0332                                       Phi phiStep_,
0333                                       Eta etaStep_) {
0334     std::vector<T> clusters;
0335     for (int phibin = 0; phibin < phiBins_; ++phibin) {  //Find eta-phibin with highest pT
0336       if (L1clusters[phibin].empty())
0337         continue;
0338
0339       // sort L1 clusters max -> min
0340       sort(L1clusters[phibin].begin(), L1clusters[phibin].end(), [](T &a, T &b) { return a.pTtot > b.pTtot; });
0341
0342       for (unsigned int imax = 0; imax < L1clusters[phibin].size(); ++imax) {
0343         if (L1clusters[phibin][imax].used)
0344           continue;
0345         Pt pt_current = L1clusters[phibin][imax].pTtot;  //current cluster (pt0)
0346         Pt pt_next = 0;                                  // next phi bin (pt1)
0347         Pt pt_next2 = 0;                                 // next to next phi bin2 (pt2)
0348         int trk1 = 0;
0349         int trk2 = 0;
0350         int tdtrk1 = 0;
0351         int tdtrk2 = 0;
0352         std::vector<unsigned int> trkidx1;
0353         std::vector<unsigned int> trkidx2;
0354         clusters.push_back(L1clusters[phibin][imax]);
0355
0356         L1clusters[phibin][imax].used = true;
0357
0358         // if we are in the last phi bin, dont check phi+1 phi+2
0359         if (phibin == phiBins_ - 1)
0360           continue;
0361
0362         std::vector<unsigned int> used_already;  //keep phi+1 clusters that have been used
0363         for (unsigned int icluster = 0; icluster < L1clusters[phibin + 1].size(); ++icluster) {
0364           if (L1clusters[phibin + 1][icluster].used)
0365             continue;
0366
0367           if (((L1clusters[phibin + 1][icluster].eta - L1clusters[phibin][imax].eta) > (3 * etaStep_) / 2) ||
0368               ((L1clusters[phibin + 1][icluster].eta - L1clusters[phibin][imax].eta) < -(3 * etaStep_) / 2))
0369             continue;
0370
0371           pt_next += L1clusters[phibin + 1][icluster].pTtot;
0372           trk1 += L1clusters[phibin + 1][icluster].ntracks;
0373           tdtrk1 += L1clusters[phibin + 1][icluster].nxtracks;
0374
0375           for (unsigned int itrk = 0; itrk < L1clusters[phibin + 1][icluster].trackidx.size(); itrk++)
0376             trkidx1.push_back(L1clusters[phibin + 1][icluster].trackidx[itrk]);
0377           used_already.push_back(icluster);
0378         }
0379
0380         if (pt_next < pt_current) {  // if pt1<pt1, merge both clusters
0381           Fill_L2Cluster<T, Pt>(clusters[clusters.size() - 1], pt_next, trk1, tdtrk1, trkidx1);
0382           for (unsigned int iused : used_already)
0383             L1clusters[phibin + 1][iused].used = true;
0384           continue;
0385         }
0386         // if phi = next to last bin there is no "next to next"
0387         if (phibin == phiBins_ - 2) {
0388           Fill_L2Cluster<T, Pt>(clusters[clusters.size() - 1], pt_next, trk1, tdtrk1, trkidx1);
0389           clusters[clusters.size() - 1].phi = L1clusters[phibin + 1][used_already[0]].phi;
0390           for (unsigned int iused : used_already)
0391             L1clusters[phibin + 1][iused].used = true;
0392           continue;
0393         }
0394         std::vector<int> used_already2;  //keep used clusters in phi+2
0395         for (unsigned int icluster = 0; icluster < L1clusters[phibin + 2].size(); ++icluster) {
0396           if (L1clusters[phibin + 2][icluster].used)
0397             continue;
0398           if (((L1clusters[phibin + 2][icluster].eta - L1clusters[phibin][imax].eta) > (3 * etaStep_) / 2) ||
0399               ((L1clusters[phibin + 2][icluster].eta - L1clusters[phibin][imax].eta) < -(3 * etaStep_) / 2))
0400             continue;
0401           pt_next2 += L1clusters[phibin + 2][icluster].pTtot;
0402           trk2 += L1clusters[phibin + 2][icluster].ntracks;
0403           tdtrk2 += L1clusters[phibin + 2][icluster].nxtracks;
0404
0405           for (unsigned int itrk = 0; itrk < L1clusters[phibin + 2][icluster].trackidx.size(); itrk++)
0406             trkidx2.push_back(L1clusters[phibin + 2][icluster].trackidx[itrk]);
0407           used_already2.push_back(icluster);
0408         }
0409         if (pt_next2 < pt_next) {
0410           std::vector<unsigned int> trkidx_both;
0411           trkidx_both.reserve(trkidx1.size() + trkidx2.size());
0412           trkidx_both.insert(trkidx_both.end(), trkidx1.begin(), trkidx1.end());
0413           trkidx_both.insert(trkidx_both.end(), trkidx2.begin(), trkidx2.end());
0414           Fill_L2Cluster<T, Pt>(
0415               clusters[clusters.size() - 1], pt_next + pt_next2, trk1 + trk2, tdtrk1 + tdtrk2, trkidx_both);
0416           clusters[clusters.size() - 1].phi = L1clusters[phibin + 1][used_already[0]].phi;
0417           for (unsigned int iused : used_already)
0418             L1clusters[phibin + 1][iused].used = true;
0419           for (unsigned int iused : used_already2)
0420             L1clusters[phibin + 2][iused].used = true;
0421         }
0422       }  // for each L1 cluster
0423     }  // for each phibin
0424
0425     int nclust = clusters.size();
0426
0427     // merge close-by clusters
0428     for (int m = 0; m < nclust - 1; ++m) {
0429       for (int n = m + 1; n < nclust; ++n) {
0430         if (clusters[n].eta != clusters[m].eta)
0431           continue;
0432         if ((DPhi(clusters[n].phi, clusters[m].phi) > (3 * phiStep_) / 2) ||
0433             (DPhi(clusters[n].phi, clusters[m].phi) < -(3 * phiStep_) / 2))
0434           continue;
0435         if (clusters[n].pTtot > clusters[m].pTtot)
0436           clusters[m].phi = clusters[n].phi;
0437         clusters[m].pTtot += clusters[n].pTtot;
0438         clusters[m].ntracks += clusters[n].ntracks;
0439         clusters[m].nxtracks += clusters[n].nxtracks;
0440         for (unsigned int itrk = 0; itrk < clusters[n].trackidx.size(); itrk++)
0441           clusters[m].trackidx.push_back(clusters[n].trackidx[itrk]);
0442         for (int m1 = n; m1 < nclust - 1; ++m1)
0443           clusters[m1] = clusters[m1 + 1];
0444         clusters.erase(clusters.begin() + nclust);
0445
0446         nclust--;
0447       }  // end of n-loop
0448     }  // end of m-loop
0449     return clusters;
0450   }
0451 }  // namespace l1ttrackjet
0452 #endif