L1TGEM/src/ME0StubAlgoPatUnit.cc

0001 #include "L1Trigger/L1TGEM/interface/ME0StubAlgoPatUnit.h"
0002
0003 using namespace l1t::me0;
0004
0005 std::vector<uint64_t> l1t::me0::maskLayerData(const std::vector<uint64_t>& data, const Mask& mask) {
0006   std::vector<uint64_t> out;
0007   out.reserve(static_cast<int>(data.size()));
0008   for (int i = 0; i < static_cast<int>(data.size()); ++i) {
0009     out.push_back(data[i] & mask.mask[i]);
0010   }
0011   return out;
0012 }
0013
0014 std::pair<std::vector<double>, double> l1t::me0::calculateCentroids(
0015     const std::vector<uint64_t>& maskedData, const std::vector<std::vector<int>>& partitionBxData) {
0016   std::vector<double> centroids;
0017   std::vector<int> bxs;
0018   for (int ly = 0; ly < static_cast<int>(maskedData.size()); ++ly) {
0019     auto data = maskedData[ly];
0020     auto bxData = partitionBxData[ly];
0021     const auto temp = findCentroid(data);
0022     double curCentroid = temp.first;
0023     std::vector<int> hitsIndices = temp.second;
0024     centroids.push_back(curCentroid);
0025
0026     for (int hitIdx : hitsIndices) {
0027       bxs.push_back(bxData[hitIdx - 1]);
0028     }
0029   }
0030   if (static_cast<int>(bxs.size()) == 0) {
0031     return {centroids, -9999};
0032   }
0033   double bxSum = std::accumulate(bxs.begin(), bxs.end(), 0.0);
0034   double count = bxs.size();
0035   return {centroids, bxSum / count};
0036 }
0037
0038 int l1t::me0::calculateHitCount(const std::vector<uint64_t>& maskedData, bool light) {
0039   int totHitCount = 0;
0040   if (light) {
0041     for (int ly : {0, 5}) {
0042       int hitLy = countOnes(maskedData[ly]);
0043       totHitCount += (hitLy < 7) ? hitLy : 7;
0044     }
0045   } else {
0046     for (uint64_t d : maskedData) {
0047       totHitCount += countOnes(d);
0048     }
0049   }
0050   return totHitCount;
0051 }
0052
0053 int l1t::me0::calculateLayerCount(const std::vector<uint64_t>& maskedData) {
0054   int lyCount = 0;
0055   bool notZero;
0056   for (uint64_t d : maskedData) {
0057     notZero = (d != 0);
0058     lyCount += static_cast<int>(notZero);
0059   }
0060   return lyCount;
0061 }
0062
0063 std::vector<int> l1t::me0::calculateClusterSize(const std::vector<uint64_t>& data) {
0064   std::vector<int> clusterSizePerLayer;
0065   clusterSizePerLayer.reserve(data.size());
0066   for (uint64_t x : data) {
0067     clusterSizePerLayer.push_back(maxClusterSize(x));
0068   }
0069   return clusterSizePerLayer;
0070 }
0071
0072 std::vector<int> l1t::me0::calculateHits(const std::vector<uint64_t>& data) {
0073   std::vector<int> nHitsPerLayer;
0074   nHitsPerLayer.reserve(data.size());
0075   for (uint64_t x : data) {
0076     nHitsPerLayer.push_back(countOnes(x));
0077   }
0078   return nHitsPerLayer;
0079 }
0080
0081 ME0StubPrimitive l1t::me0::patUnit(const std::vector<uint64_t>& data,
0082                                    const std::vector<std::vector<int>>& bxData,
0083                                    int strip,
0084                                    int partition,
0085                                    const std::vector<int>& lyThreshPatid,
0086                                    const std::vector<int>& lyThreshEta,
0087                                    int inputMaxSpan,
0088                                    bool skipCentroids,
0089                                    int numOr,
0090                                    bool lightHitCount,
0091                                    bool verbose) {
0092   // construct the dynamic_patlist (we do not use default PATLIST anymore)
0093   // for robustness concern, other codes might use PATLIST, so we kept the default PATLIST in subfunc
0094   // however, this could cause inconsistent issue, becareful! OR find a way to modify PATLIST
0095
0096   /*
0097     takes in sample data for each layer and returns best segment
0098
0099     processing pipeline is
0100
0101     (1) take in 6 layers of raw data
0102     (2) for the X (~16) patterns available, AND together the raw data with the respective pattern masks
0103     (3) count the # of hits in each pattern
0104     (4) calculate the centroids for each pattern
0105     (5) process segments
0106     (6) choose the max of all patterns
0107     (7) apply a layer threshold
0108     */
0109
0110   // (2)
0111   // and the layer data with the respective layer mask to
0112   // determine how many hits are in each layer
0113   // this yields a map object that can be iterated over to get,
0114   //    for each of the N patterns, the masked []*6 layer data
0115   std::vector<std::vector<uint64_t>> maskedData;
0116   std::vector<int> pids;
0117   for (const Mask& M : kLayerMask) {
0118     maskedData.push_back(maskLayerData(data, M));
0119     pids.push_back(M.id);
0120   }
0121
0122   // (3) count # of hits & process centroids
0123   std::vector<int> hcs;
0124   std::vector<int> lcs;
0125   std::vector<std::vector<double>> centroids;
0126   std::vector<double> bxs;
0127   for (const std::vector<uint64_t>& x : maskedData) {
0128     hcs.push_back(calculateHitCount(x, lightHitCount));
0129     lcs.push_back(calculateLayerCount(x));
0130     if (skipCentroids) {
0131       centroids.push_back({0, 0, 0, 0, 0, 0});
0132       bxs.push_back(-9999);
0133     } else {
0134       auto temp = calculateCentroids(x, bxData);
0135       std::vector<double> curPatternCentroids = temp.first;
0136       int curPatternBx = temp.second;
0137       centroids.push_back(curPatternCentroids);
0138       bxs.push_back(curPatternBx);
0139     }
0140   }
0141
0142   // (4) process segments & choose the max of all patterns
0143   ME0StubPrimitive best{0, 0, 0, strip, partition};
0144   for (int i = 0; i < static_cast<int>(hcs.size()); ++i) {
0145     ME0StubPrimitive seg{lcs[i], hcs[i], pids[i], strip, partition, bxs[i]};
0146     seg.updateQuality();
0147     if (best.quality() < seg.quality()) {
0148       best = seg;
0149       best.setCentroids(centroids[i]);
0150       best.updateQuality();
0151     }
0152   }
0153
0154   // (5) apply a layer threshold
0155   int lyThreshFinal;
0156   if (lyThreshPatid[best.patternId() - 1] > lyThreshEta[partition]) {
0157     lyThreshFinal = lyThreshPatid[best.patternId() - 1];
0158   } else {
0159     lyThreshFinal = lyThreshEta[partition];
0160   }
0161
0162   if (best.layerCount() < lyThreshFinal) {
0163     best.reset();
0164   }
0165
0166   // (6) remove very wide segments
0167   if (best.patternId() <= 10) {
0168     best.reset();
0169   }
0170
0171   // (7) remove segments with large clusters for wide segments - ONLY NEEDED FOR PU200 - NOT USED AT THE MOMENT
0172   std::vector<int> clusterSizeMaxLimits = {3, 6, 9, 12, 15};
0173   std::vector<int> nHitsMaxLimits = {3, 6, 9, 12, 15};
0174   std::vector<int> clusterSizeCounts = calculateClusterSize(data);
0175   std::vector<int> nHitsCounts = calculateHits(data);
0176   std::vector<int> nLayersLargeClusters = {0, 0, 0, 0, 0};
0177   std::vector<int> nLayersLargeHits = {0, 0, 0, 0, 0};
0178   for (int i = 0; i < static_cast<int>(clusterSizeCounts.size()); ++i) {
0179     int threshold = clusterSizeMaxLimits[i];
0180     for (int l : clusterSizeCounts) {
0181       if (l > threshold) {
0182         nLayersLargeClusters[i]++;
0183       }
0184     }
0185   }
0186   for (int i = 0; i < static_cast<int>(nHitsMaxLimits.size()); ++i) {
0187     int threshold = nHitsMaxLimits[i];
0188     for (int l : nHitsCounts) {
0189       if (l > threshold) {
0190         nLayersLargeHits[i]++;
0191       }
0192     }
0193   }
0194
0195   best.setMaxClusterSize(*std::max_element(clusterSizeCounts.begin(), clusterSizeCounts.end()));
0196   best.setMaxNoise(*std::max_element(nHitsCounts.begin(), nHitsCounts.end()));
0197
0198   best.setHitCount(0);
0199   best.updateQuality();
0200
0201   return best;
0202 }