Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-05-12 02:41:22 

0001 #include "L1Trigger/L1THGCal/interface/backend/HGCalHistoSeedingImpl.h"
0002 #include "L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h"
0003 #include "DataFormats/Math/interface/deltaR.h"
0004 #include "DataFormats/Math/interface/deltaPhi.h"
0005 #include <numeric>
0006 
0007 HGCalHistoSeedingImpl::HGCalHistoSeedingImpl(const edm::ParameterSet& conf)
0008     : seedingAlgoType_(conf.getParameter<std::string>("type_histoalgo")),
0009       nBins1_(conf.getParameter<unsigned>("nBins_X1_histo_multicluster")),
0010       nBins2_(conf.getParameter<unsigned>("nBins_X2_histo_multicluster")),
0011       binsSumsHisto_(conf.getParameter<std::vector<unsigned>>("binSumsHisto")),
0012       histoThreshold_(conf.getParameter<double>("threshold_histo_multicluster")),
0013       neighbour_weights_(conf.getParameter<std::vector<double>>("neighbour_weights")),
0014       smoothing_ecal_(conf.getParameter<std::vector<double>>("seed_smoothing_ecal")),
0015       smoothing_hcal_(conf.getParameter<std::vector<double>>("seed_smoothing_hcal")),
0016       kROverZMin_(conf.getParameter<double>("kROverZMin")),
0017       kROverZMax_(conf.getParameter<double>("kROverZMax")) {
0018   if (seedingAlgoType_ == "HistoMaxC3d") {
0019     seedingType_ = HistoMaxC3d;
0020   } else if (seedingAlgoType_ == "HistoSecondaryMaxC3d") {
0021     seedingType_ = HistoSecondaryMaxC3d;
0022   } else if (seedingAlgoType_ == "HistoThresholdC3d") {
0023     seedingType_ = HistoThresholdC3d;
0024   } else if (seedingAlgoType_ == "HistoInterpolatedMaxC3d") {
0025     seedingType_ = HistoInterpolatedMaxC3d;
0026   } else {
0027     throw cms::Exception("HGCTriggerParameterError") << "Unknown Multiclustering type '" << seedingAlgoType_;
0028   }
0029 
0030   if (conf.getParameter<std::string>("seed_position") == "BinCentre") {
0031     seedingPosition_ = BinCentre;
0032   } else if (conf.getParameter<std::string>("seed_position") == "TCWeighted") {
0033     seedingPosition_ = TCWeighted;
0034   } else {
0035     throw cms::Exception("HGCTriggerParameterError")
0036         << "Unknown Seed Position option '" << conf.getParameter<std::string>("seed_position");
0037   }
0038   if (conf.getParameter<std::string>("seeding_space") == "RPhi") {
0039     seedingSpace_ = RPhi;
0040     navigator_ = Navigator(nBins1_, Navigator::AxisType::Bounded, nBins2_, Navigator::AxisType::Circular);
0041   } else if (conf.getParameter<std::string>("seeding_space") == "XY") {
0042     seedingSpace_ = XY;
0043     navigator_ = Navigator(nBins1_, Navigator::AxisType::Bounded, nBins2_, Navigator::AxisType::Bounded);
0044   } else {
0045     throw cms::Exception("HGCTriggerParameterError")
0046         << "Unknown seeding space  '" << conf.getParameter<std::string>("seeding_space");
0047   }
0048 
0049   edm::LogInfo("HGCalMulticlusterParameters")
0050       << "\nMulticluster number of X1-bins for the histo algorithm: " << nBins1_
0051       << "\nMulticluster number of X2-bins for the histo algorithm: " << nBins2_
0052       << "\nMulticluster MIPT threshold for histo threshold algorithm: " << histoThreshold_
0053       << "\nMulticluster type of multiclustering algortihm: " << seedingAlgoType_;
0054 
0055   if (seedingAlgoType_.find("Histo") != std::string::npos && seedingSpace_ == RPhi &&
0056       nBins1_ != binsSumsHisto_.size()) {
0057     throw cms::Exception("Inconsistent bin size")
0058         << "Inconsistent nBins_X1_histo_multicluster ( " << nBins1_ << " ) and binSumsHisto ( " << binsSumsHisto_.size()
0059         << " ) size in HGCalMulticlustering\n";
0060   }
0061 
0062   if (neighbour_weights_.size() != neighbour_weights_size_) {
0063     throw cms::Exception("Inconsistent vector size")
0064         << "Inconsistent size of neighbour weights vector in HGCalMulticlustering ( " << neighbour_weights_.size()
0065         << " ). Should be " << neighbour_weights_size_ << "\n";
0066   }
0067 }
0068 
0069 HGCalHistoSeedingImpl::Histogram HGCalHistoSeedingImpl::fillHistoClusters(
0070     const std::vector<edm::Ptr<l1t::HGCalCluster>>& clustersPtrs) {
0071   Histogram histoClusters(nBins1_, nBins2_);
0072   std::array<double, 4> bounds = boundaries();
0073   double minx1 = std::get<0>(bounds);
0074   double maxx1 = std::get<1>(bounds);
0075   double minx2 = std::get<2>(bounds);
0076   double maxx2 = std::get<3>(bounds);
0077 
0078   for (auto& clu : clustersPtrs) {
0079     float x1 = 0., x2 = 0;
0080     switch (seedingSpace_) {
0081       case RPhi:
0082         x1 = sqrt(pow(clu->centreProj().x(), 2) + pow(clu->centreProj().y(), 2));
0083         x2 = reco::reduceRange(clu->phi());
0084         break;
0085       case XY:
0086         x1 = clu->centreProj().x();
0087         x2 = clu->centreProj().y();
0088         break;
0089     };
0090     if (x1 < minx1 || x1 >= maxx1) {
0091       throw cms::Exception("OutOfBound") << "TC X1 = " << x1 << " out of the seeding histogram bounds " << minx1
0092                                          << " - " << maxx1;
0093     }
0094     if (x2 < minx2 || x2 >= maxx2) {
0095       throw cms::Exception("OutOfBound") << "TC X2 = " << x2 << " out of the seeding histogram bounds " << minx2
0096                                          << " - " << maxx2;
0097     }
0098     unsigned bin1 = unsigned((x1 - minx1) * nBins1_ / (maxx1 - minx1));
0099     unsigned bin2 = unsigned((x2 - minx2) * nBins2_ / (maxx2 - minx2));
0100 
0101     auto& bin = histoClusters.at(triggerTools_.zside(clu->detId()), bin1, bin2);
0102     bin.values[Bin::Content::Sum] += clu->mipPt();
0103     if (triggerTools_.isEm(clu->detId())) {
0104       bin.values[Bin::Content::Ecal] += clu->mipPt();
0105     } else {
0106       bin.values[Bin::Content::Hcal] += clu->mipPt();
0107     }
0108     bin.weighted_x += (clu->centreProj().x()) * clu->mipPt();
0109     bin.weighted_y += (clu->centreProj().y()) * clu->mipPt();
0110   }
0111 
0112   for (auto& bin : histoClusters) {
0113     bin.weighted_x /= bin.values[Bin::Content::Sum];
0114     bin.weighted_y /= bin.values[Bin::Content::Sum];
0115   }
0116 
0117   return histoClusters;
0118 }
0119 
0120 HGCalHistoSeedingImpl::Histogram HGCalHistoSeedingImpl::fillSmoothHistoClusters(const Histogram& histoClusters,
0121                                                                                 const vector<double>& kernel,
0122                                                                                 Bin::Content binContent) {
0123   Histogram histoSmooth(histoClusters);
0124 
0125   unsigned kernel_size = std::sqrt(kernel.size());
0126   if (kernel_size * kernel_size != kernel.size()) {
0127     throw cms::Exception("HGCTriggerParameterError") << "Only square kernels can be used.";
0128   }
0129   if (kernel_size % 2 != 1) {
0130     throw cms::Exception("HGCTriggerParameterError") << "The kernel size must be an odd value.";
0131   }
0132   int shift_max = (kernel_size - 1) / 2;
0133   double normalization = std::accumulate(kernel.begin(), kernel.end(), 0.);
0134   for (int z_side : {-1, 1}) {
0135     for (unsigned x1 = 0; x1 < nBins1_; x1++) {
0136       for (unsigned x2 = 0; x2 < nBins2_; x2++) {
0137         const auto& bin_orig = histoClusters.at(z_side, x1, x2);
0138         double smooth = 0.;
0139         navigator_.setHome(x1, x2);
0140         for (int x1_shift = -shift_max; x1_shift <= shift_max; x1_shift++) {
0141           int index1 = x1_shift + shift_max;
0142           for (int x2_shift = -shift_max; x2_shift <= shift_max; x2_shift++) {
0143             auto shifted = navigator_.move(x1_shift, x2_shift);
0144             int index2 = x2_shift + shift_max;
0145             double kernel_value = kernel.at(index1 * kernel_size + index2);
0146             bool out = shifted[0] == -1 || shifted[1] == -1;
0147             double content = (out ? 0. : histoClusters.at(z_side, shifted[0], shifted[1]).values[binContent]);
0148             smooth += content * kernel_value;
0149           }
0150         }
0151         auto& bin = histoSmooth.at(z_side, x1, x2);
0152         bin.values[binContent] = smooth / normalization;
0153         bin.weighted_x = bin_orig.weighted_x;
0154         bin.weighted_y = bin_orig.weighted_y;
0155       }
0156     }
0157   }
0158 
0159   return histoSmooth;
0160 }
0161 
0162 HGCalHistoSeedingImpl::Histogram HGCalHistoSeedingImpl::fillSmoothPhiHistoClusters(const Histogram& histoClusters,
0163                                                                                    const vector<unsigned>& binSums) {
0164   Histogram histoSumPhiClusters(nBins1_, nBins2_);
0165 
0166   for (int z_side : {-1, 1}) {
0167     for (unsigned bin1 = 0; bin1 < nBins1_; bin1++) {
0168       int nBinsSide = (binSums[bin1] - 1) / 2;
0169       float R1 = kROverZMin_ + bin1 * (kROverZMax_ - kROverZMin_) / nBins1_;
0170       float R2 = R1 + ((kROverZMax_ - kROverZMin_) / nBins1_);
0171       double area =
0172           ((M_PI * (pow(R2, 2) - pow(R1, 2))) / nBins2_) *
0173           (1 +
0174            2.0 *
0175                (1 -
0176                 pow(0.5,
0177                     nBinsSide)));  // Takes into account different area of bins in different R-rings + sum of quadratic weights used
0178 
0179       for (unsigned bin2 = 0; bin2 < nBins2_; bin2++) {
0180         const auto& bin_orig = histoClusters.at(z_side, bin1, bin2);
0181         float content = bin_orig.values[Bin::Content::Sum];
0182 
0183         for (int bin22 = 1; bin22 <= nBinsSide; bin22++) {
0184           int binToSumLeft = bin2 - bin22;
0185           if (binToSumLeft < 0)
0186             binToSumLeft += nBins2_;
0187           unsigned binToSumRight = bin2 + bin22;
0188           if (binToSumRight >= nBins2_)
0189             binToSumRight -= nBins2_;
0190 
0191           content += histoClusters.at(z_side, bin1, binToSumLeft).values[Bin::Content::Sum] /
0192                      pow(2, bin22);  // quadratic kernel
0193 
0194           content += histoClusters.at(z_side, bin1, binToSumRight).values[Bin::Content::Sum] /
0195                      pow(2, bin22);  // quadratic kernel
0196         }
0197 
0198         auto& bin = histoSumPhiClusters.at(z_side, bin1, bin2);
0199         bin.values[Bin::Content::Sum] = (content / area) * area_per_triggercell_;
0200         bin.weighted_x = bin_orig.weighted_x;
0201         bin.weighted_y = bin_orig.weighted_y;
0202       }
0203     }
0204   }
0205 
0206   return histoSumPhiClusters;
0207 }
0208 
0209 HGCalHistoSeedingImpl::Histogram HGCalHistoSeedingImpl::fillSmoothRPhiHistoClusters(const Histogram& histoClusters) {
0210   Histogram histoSumRPhiClusters(nBins1_, nBins2_);
0211 
0212   for (int z_side : {-1, 1}) {
0213     for (unsigned bin1 = 0; bin1 < nBins1_; bin1++) {
0214       float weight =
0215           (bin1 == 0 || bin1 == nBins1_ - 1) ? 1.5 : 2.;  //Take into account edges with only one side up or down
0216 
0217       for (unsigned bin2 = 0; bin2 < nBins2_; bin2++) {
0218         const auto& bin_orig = histoClusters.at(z_side, bin1, bin2);
0219         float content = bin_orig.values[Bin::Content::Sum];
0220         float contentDown = bin1 > 0 ? histoClusters.at(z_side, bin1 - 1, bin2).values[Bin::Content::Sum] : 0;
0221         float contentUp = bin1 < (nBins1_ - 1) ? histoClusters.at(z_side, bin1 + 1, bin2).values[Bin::Content::Sum] : 0;
0222 
0223         auto& bin = histoSumRPhiClusters.at(z_side, bin1, bin2);
0224         bin.values[Bin::Content::Sum] = (content + 0.5 * contentDown + 0.5 * contentUp) / weight;
0225         bin.weighted_x = bin_orig.weighted_x;
0226         bin.weighted_y = bin_orig.weighted_y;
0227       }
0228     }
0229   }
0230 
0231   return histoSumRPhiClusters;
0232 }
0233 
0234 void HGCalHistoSeedingImpl::setSeedEnergyAndPosition(std::vector<std::pair<GlobalPoint, double>>& seedPositionsEnergy,
0235                                                      int z_side,
0236                                                      unsigned bin1,
0237                                                      unsigned bin2,
0238                                                      const Bin& histBin) {
0239   float x_seed = 0;
0240   float y_seed = 0;
0241   std::array<double, 4> bounds = boundaries();
0242   double minx1 = std::get<0>(bounds);
0243   double maxx1 = std::get<1>(bounds);
0244   double minx2 = std::get<2>(bounds);
0245   double maxx2 = std::get<3>(bounds);
0246 
0247   if (seedingPosition_ == BinCentre) {
0248     float x1_seed = minx1 + (bin1 + 0.5) * (maxx1 - minx1) / nBins1_;
0249     float x2_seed = minx2 + (bin2 + 0.5) * (maxx2 - minx2) / nBins2_;
0250     switch (seedingSpace_) {
0251       case RPhi:
0252         x_seed = x1_seed * cos(x2_seed);
0253         y_seed = x1_seed * sin(x2_seed);
0254         break;
0255       case XY:
0256         x_seed = x1_seed;
0257         y_seed = x2_seed;
0258     };
0259   } else if (seedingPosition_ == TCWeighted) {
0260     x_seed = histBin.weighted_x;
0261     y_seed = histBin.weighted_y;
0262   }
0263 
0264   seedPositionsEnergy.emplace_back(GlobalPoint(x_seed, y_seed, z_side), histBin.values[Bin::Content::Sum]);
0265 }
0266 
0267 std::vector<std::pair<GlobalPoint, double>> HGCalHistoSeedingImpl::computeMaxSeeds(const Histogram& histoClusters) {
0268   std::vector<std::pair<GlobalPoint, double>> seedPositionsEnergy;
0269 
0270   for (int z_side : {-1, 1}) {
0271     for (unsigned bin1 = 0; bin1 < nBins1_; bin1++) {
0272       for (unsigned bin2 = 0; bin2 < nBins2_; bin2++) {
0273         float MIPT_seed = histoClusters.at(z_side, bin1, bin2).values[Bin::Content::Sum];
0274         bool isMax = MIPT_seed > histoThreshold_;
0275         if (!isMax)
0276           continue;
0277 
0278         navigator_.setHome(bin1, bin2);
0279         auto pos_N = navigator_.move(1, 0);
0280         auto pos_S = navigator_.move(-1, 0);
0281         auto pos_W = navigator_.move(0, -1);
0282         auto pos_E = navigator_.move(0, 1);
0283         auto pos_NW = navigator_.move(1, -1);
0284         auto pos_NE = navigator_.move(1, 1);
0285         auto pos_SW = navigator_.move(-1, -1);
0286         auto pos_SE = navigator_.move(-1, 1);
0287 
0288         float MIPT_N = (pos_N[0] != -1 && pos_N[1] != -1)
0289                            ? histoClusters.at(z_side, pos_N[0], pos_N[1]).values[Bin::Content::Sum]
0290                            : 0;
0291         float MIPT_S = (pos_S[0] != -1 && pos_S[1] != -1)
0292                            ? histoClusters.at(z_side, pos_S[0], pos_S[1]).values[Bin::Content::Sum]
0293                            : 0;
0294         float MIPT_W = (pos_W[0] != -1 && pos_W[1] != -1)
0295                            ? histoClusters.at(z_side, pos_W[0], pos_W[1]).values[Bin::Content::Sum]
0296                            : 0;
0297         float MIPT_E = (pos_E[0] != -1 && pos_E[1] != -1)
0298                            ? histoClusters.at(z_side, pos_E[0], pos_E[1]).values[Bin::Content::Sum]
0299                            : 0;
0300         float MIPT_NW = (pos_NW[0] != -1 && pos_NW[1] != -1)
0301                             ? histoClusters.at(z_side, pos_NW[0], pos_NW[1]).values[Bin::Content::Sum]
0302                             : 0;
0303         float MIPT_NE = (pos_NE[0] != -1 && pos_NE[1] != -1)
0304                             ? histoClusters.at(z_side, pos_NE[0], pos_NE[1]).values[Bin::Content::Sum]
0305                             : 0;
0306         float MIPT_SW = (pos_SW[0] != -1 && pos_SW[1] != -1)
0307                             ? histoClusters.at(z_side, pos_SW[0], pos_SW[1]).values[Bin::Content::Sum]
0308                             : 0;
0309         float MIPT_SE = (pos_SE[0] != -1 && pos_SE[1] != -1)
0310                             ? histoClusters.at(z_side, pos_SE[0], pos_SE[1]).values[Bin::Content::Sum]
0311                             : 0;
0312 
0313         isMax &= MIPT_seed >= MIPT_S && MIPT_seed > MIPT_N && MIPT_seed >= MIPT_E && MIPT_seed >= MIPT_SE &&
0314                  MIPT_seed >= MIPT_NE && MIPT_seed > MIPT_W && MIPT_seed > MIPT_SW && MIPT_seed > MIPT_NW;
0315 
0316         if (isMax) {
0317           setSeedEnergyAndPosition(seedPositionsEnergy, z_side, bin1, bin2, histoClusters.at(z_side, bin1, bin2));
0318         }
0319       }
0320     }
0321   }
0322 
0323   return seedPositionsEnergy;
0324 }
0325 
0326 std::vector<std::pair<GlobalPoint, double>> HGCalHistoSeedingImpl::computeInterpolatedMaxSeeds(
0327     const Histogram& histoClusters) {
0328   std::vector<std::pair<GlobalPoint, double>> seedPositionsEnergy;
0329 
0330   for (int z_side : {-1, 1}) {
0331     for (unsigned bin1 = 0; bin1 < nBins1_; bin1++) {
0332       for (unsigned bin2 = 0; bin2 < nBins2_; bin2++) {
0333         float MIPT_seed = histoClusters.at(z_side, bin1, bin2).values[Bin::Content::Sum];
0334 
0335         navigator_.setHome(bin1, bin2);
0336         auto pos_N = navigator_.move(1, 0);
0337         auto pos_S = navigator_.move(-1, 0);
0338         auto pos_W = navigator_.move(0, -1);
0339         auto pos_E = navigator_.move(0, 1);
0340         auto pos_NW = navigator_.move(1, -1);
0341         auto pos_NE = navigator_.move(1, 1);
0342         auto pos_SW = navigator_.move(-1, -1);
0343         auto pos_SE = navigator_.move(-1, 1);
0344 
0345         float MIPT_N = (pos_N[0] != -1 && pos_N[1] != -1)
0346                            ? histoClusters.at(z_side, pos_N[0], pos_N[1]).values[Bin::Content::Sum]
0347                            : 0;
0348         float MIPT_S = (pos_S[0] != -1 && pos_S[1] != -1)
0349                            ? histoClusters.at(z_side, pos_S[0], pos_S[1]).values[Bin::Content::Sum]
0350                            : 0;
0351         float MIPT_W = (pos_W[0] != -1 && pos_W[1] != -1)
0352                            ? histoClusters.at(z_side, pos_W[0], pos_W[1]).values[Bin::Content::Sum]
0353                            : 0;
0354         float MIPT_E = (pos_E[0] != -1 && pos_E[1] != -1)
0355                            ? histoClusters.at(z_side, pos_E[0], pos_E[1]).values[Bin::Content::Sum]
0356                            : 0;
0357         float MIPT_NW = (pos_NW[0] != -1 && pos_NW[1] != -1)
0358                             ? histoClusters.at(z_side, pos_NW[0], pos_NW[1]).values[Bin::Content::Sum]
0359                             : 0;
0360         float MIPT_NE = (pos_NE[0] != -1 && pos_NE[1] != -1)
0361                             ? histoClusters.at(z_side, pos_NE[0], pos_NE[1]).values[Bin::Content::Sum]
0362                             : 0;
0363         float MIPT_SW = (pos_SW[0] != -1 && pos_SW[1] != -1)
0364                             ? histoClusters.at(z_side, pos_SW[0], pos_SW[1]).values[Bin::Content::Sum]
0365                             : 0;
0366         float MIPT_SE = (pos_SE[0] != -1 && pos_SE[1] != -1)
0367                             ? histoClusters.at(z_side, pos_SE[0], pos_SE[1]).values[Bin::Content::Sum]
0368                             : 0;
0369 
0370         float MIPT_pred = neighbour_weights_.at(0) * MIPT_NW + neighbour_weights_.at(1) * MIPT_N +
0371                           neighbour_weights_.at(2) * MIPT_NE + neighbour_weights_.at(3) * MIPT_W +
0372                           neighbour_weights_.at(5) * MIPT_E + neighbour_weights_.at(6) * MIPT_SW +
0373                           neighbour_weights_.at(7) * MIPT_S + neighbour_weights_.at(8) * MIPT_SE;
0374 
0375         bool isMax = MIPT_seed >= (MIPT_pred + histoThreshold_);
0376 
0377         if (isMax) {
0378           setSeedEnergyAndPosition(seedPositionsEnergy, z_side, bin1, bin2, histoClusters.at(z_side, bin1, bin2));
0379         }
0380       }
0381     }
0382   }
0383 
0384   return seedPositionsEnergy;
0385 }
0386 
0387 std::vector<std::pair<GlobalPoint, double>> HGCalHistoSeedingImpl::computeThresholdSeeds(
0388     const Histogram& histoClusters) {
0389   std::vector<std::pair<GlobalPoint, double>> seedPositionsEnergy;
0390 
0391   for (int z_side : {-1, 1}) {
0392     for (unsigned bin1 = 0; bin1 < nBins1_; bin1++) {
0393       for (unsigned bin2 = 0; bin2 < nBins2_; bin2++) {
0394         float MIPT_seed = histoClusters.at(z_side, bin1, bin2).values[Bin::Content::Sum];
0395         bool isSeed = MIPT_seed > histoThreshold_;
0396 
0397         if (isSeed) {
0398           setSeedEnergyAndPosition(seedPositionsEnergy, z_side, bin1, bin2, histoClusters.at(z_side, bin1, bin2));
0399         }
0400       }
0401     }
0402   }
0403 
0404   return seedPositionsEnergy;
0405 }
0406 
0407 std::vector<std::pair<GlobalPoint, double>> HGCalHistoSeedingImpl::computeSecondaryMaxSeeds(
0408     const Histogram& histoClusters) {
0409   std::vector<std::pair<GlobalPoint, double>> seedPositionsEnergy;
0410 
0411   HistogramT<uint8_t> primarySeedPositions(nBins1_, nBins2_);
0412   HistogramT<uint8_t> vetoPositions(nBins1_, nBins2_);
0413 
0414   //Search for primary seeds
0415   for (int z_side : {-1, 1}) {
0416     for (unsigned bin1 = 0; bin1 < nBins1_; bin1++) {
0417       for (unsigned bin2 = 0; bin2 < nBins2_; bin2++) {
0418         float MIPT_seed = histoClusters.at(z_side, bin1, bin2).values[Bin::Content::Sum];
0419         bool isMax = MIPT_seed > histoThreshold_;
0420 
0421         if (!isMax)
0422           continue;
0423 
0424         float MIPT_S = bin1 < (nBins1_ - 1) ? histoClusters.at(z_side, bin1 + 1, bin2).values[Bin::Content::Sum] : 0;
0425         float MIPT_N = bin1 > 0 ? histoClusters.at(z_side, bin1 - 1, bin2).values[Bin::Content::Sum] : 0;
0426 
0427         int binLeft = bin2 - 1;
0428         if (binLeft < 0)
0429           binLeft += nBins2_;
0430         unsigned binRight = bin2 + 1;
0431         if (binRight >= nBins2_)
0432           binRight -= nBins2_;
0433 
0434         float MIPT_W = histoClusters.at(z_side, bin1, binLeft).values[Bin::Content::Sum];
0435         float MIPT_E = histoClusters.at(z_side, bin1, binRight).values[Bin::Content::Sum];
0436         float MIPT_NW = bin1 > 0 ? histoClusters.at(z_side, bin1 - 1, binLeft).values[Bin::Content::Sum] : 0;
0437         float MIPT_NE = bin1 > 0 ? histoClusters.at(z_side, bin1 - 1, binRight).values[Bin::Content::Sum] : 0;
0438         float MIPT_SW =
0439             bin1 < (nBins1_ - 1) ? histoClusters.at(z_side, bin1 + 1, binLeft).values[Bin::Content::Sum] : 0;
0440         float MIPT_SE =
0441             bin1 < (nBins1_ - 1) ? histoClusters.at(z_side, bin1 + 1, binRight).values[Bin::Content::Sum] : 0;
0442 
0443         isMax &= MIPT_seed >= MIPT_S && MIPT_seed > MIPT_N && MIPT_seed >= MIPT_E && MIPT_seed >= MIPT_SE &&
0444                  MIPT_seed >= MIPT_NE && MIPT_seed > MIPT_W && MIPT_seed > MIPT_SW && MIPT_seed > MIPT_NW;
0445 
0446         if (isMax) {
0447           setSeedEnergyAndPosition(seedPositionsEnergy, z_side, bin1, bin2, histoClusters.at(z_side, bin1, bin2));
0448 
0449           primarySeedPositions.at(z_side, bin1, bin2) = true;
0450 
0451           vetoPositions.at(z_side, bin1, binLeft) = true;
0452           vetoPositions.at(z_side, bin1, binRight) = true;
0453           if (bin1 > 0) {
0454             vetoPositions.at(z_side, bin1 - 1, bin2) = true;
0455             vetoPositions.at(z_side, bin1 - 1, binRight) = true;
0456             vetoPositions.at(z_side, bin1 - 1, binLeft) = true;
0457           }
0458           if (bin1 < (nBins1_ - 1)) {
0459             vetoPositions.at(z_side, bin1 + 1, bin2) = true;
0460             vetoPositions.at(z_side, bin1 + 1, binRight) = true;
0461             vetoPositions.at(z_side, bin1 + 1, binLeft) = true;
0462           }
0463         }
0464       }
0465     }
0466   }
0467 
0468   //Search for secondary seeds
0469 
0470   for (int z_side : {-1, 1}) {
0471     for (unsigned bin1 = 0; bin1 < nBins1_; bin1++) {
0472       for (unsigned bin2 = 0; bin2 < nBins2_; bin2++) {
0473         //Cannot be a secondary seed if already a primary seed, or adjacent to primary seed
0474         if (primarySeedPositions.at(z_side, bin1, bin2) || vetoPositions.at(z_side, bin1, bin2))
0475           continue;
0476 
0477         float MIPT_seed = histoClusters.at(z_side, bin1, bin2).values[Bin::Content::Sum];
0478         bool isMax = MIPT_seed > histoThreshold_;
0479 
0480         float MIPT_S = bin1 < (nBins1_ - 1) ? histoClusters.at(z_side, bin1 + 1, bin2).values[Bin::Content::Sum] : 0;
0481         float MIPT_N = bin1 > 0 ? histoClusters.at(z_side, bin1 - 1, bin2).values[Bin::Content::Sum] : 0;
0482 
0483         int binLeft = bin2 - 1;
0484         if (binLeft < 0)
0485           binLeft += nBins2_;
0486         unsigned binRight = bin2 + 1;
0487         if (binRight >= nBins2_)
0488           binRight -= nBins2_;
0489 
0490         float MIPT_W = histoClusters.at(z_side, bin1, binLeft).values[Bin::Content::Sum];
0491         float MIPT_E = histoClusters.at(z_side, bin1, binRight).values[Bin::Content::Sum];
0492         float MIPT_NW = bin1 > 0 ? histoClusters.at(z_side, bin1 - 1, binLeft).values[Bin::Content::Sum] : 0;
0493         float MIPT_NE = bin1 > 0 ? histoClusters.at(z_side, bin1 - 1, binRight).values[Bin::Content::Sum] : 0;
0494         float MIPT_SW =
0495             bin1 < (nBins1_ - 1) ? histoClusters.at(z_side, bin1 + 1, binLeft).values[Bin::Content::Sum] : 0;
0496         float MIPT_SE =
0497             bin1 < (nBins1_ - 1) ? histoClusters.at(z_side, bin1 + 1, binRight).values[Bin::Content::Sum] : 0;
0498 
0499         isMax &= (((bin1 < nBins1_ - 1) && vetoPositions.at(z_side, bin1 + 1, bin2)) or MIPT_seed >= MIPT_S) &&
0500                  (((bin1 > 0) && vetoPositions.at(z_side, bin1 - 1, bin2)) or MIPT_seed > MIPT_N) &&
0501                  ((vetoPositions.at(z_side, bin1, binRight)) or MIPT_seed >= MIPT_E) &&
0502                  (((bin1 < nBins1_ - 1) && vetoPositions.at(z_side, bin1 + 1, binRight)) or MIPT_seed >= MIPT_SE) &&
0503                  (((bin1 > 0) && vetoPositions.at(z_side, bin1 - 1, binRight)) or MIPT_seed >= MIPT_NE) &&
0504                  ((vetoPositions.at(z_side, bin1, binLeft)) or MIPT_seed > MIPT_W) &&
0505                  (((bin1 < nBins1_ - 1) && vetoPositions.at(z_side, bin1 + 1, binLeft)) or MIPT_seed > MIPT_SW) &&
0506                  (((bin1 > 0) && vetoPositions.at(z_side, bin1 - 1, binLeft)) or MIPT_seed > MIPT_NW);
0507 
0508         if (isMax) {
0509           setSeedEnergyAndPosition(seedPositionsEnergy, z_side, bin1, bin2, histoClusters.at(z_side, bin1, bin2));
0510         }
0511       }
0512     }
0513   }
0514 
0515   return seedPositionsEnergy;
0516 }
0517 
0518 void HGCalHistoSeedingImpl::findHistoSeeds(const std::vector<edm::Ptr<l1t::HGCalCluster>>& clustersPtrs,
0519                                            std::vector<std::pair<GlobalPoint, double>>& seedPositionsEnergy) {
0520   /* put clusters into an r/z x phi histogram */
0521   Histogram histoCluster = fillHistoClusters(clustersPtrs);
0522 
0523   Histogram smoothHistoCluster;
0524   if (seedingSpace_ == RPhi) {
0525     /* smoothen along the phi direction + normalize each bin to same area */
0526     Histogram smoothPhiHistoCluster = fillSmoothPhiHistoClusters(histoCluster, binsSumsHisto_);
0527 
0528     /* smoothen along the r/z direction */
0529     smoothHistoCluster = fillSmoothRPhiHistoClusters(smoothPhiHistoCluster);
0530   } else if (seedingSpace_ == XY) {
0531     smoothHistoCluster = fillSmoothHistoClusters(histoCluster, smoothing_ecal_, Bin::Content::Ecal);
0532     smoothHistoCluster = fillSmoothHistoClusters(smoothHistoCluster, smoothing_hcal_, Bin::Content::Hcal);
0533     // Update sum with smoothen ECAL + HCAL
0534     for (int z_side : {-1, 1}) {
0535       for (unsigned x1 = 0; x1 < nBins1_; x1++) {
0536         for (unsigned x2 = 0; x2 < nBins2_; x2++) {
0537           auto& bin = smoothHistoCluster.at(z_side, x1, x2);
0538           bin.values[Bin::Content::Sum] = bin.values[Bin::Content::Ecal] + bin.values[Bin::Content::Hcal];
0539         }
0540       }
0541     }
0542   }
0543 
0544   /* seeds determined with local maximum criteria */
0545   if (seedingType_ == HistoMaxC3d)
0546     seedPositionsEnergy = computeMaxSeeds(smoothHistoCluster);
0547   else if (seedingType_ == HistoThresholdC3d)
0548     seedPositionsEnergy = computeThresholdSeeds(smoothHistoCluster);
0549   else if (seedingType_ == HistoInterpolatedMaxC3d)
0550     seedPositionsEnergy = computeInterpolatedMaxSeeds(smoothHistoCluster);
0551   else if (seedingType_ == HistoSecondaryMaxC3d)
0552     seedPositionsEnergy = computeSecondaryMaxSeeds(smoothHistoCluster);
0553 }
0554 
0555 std::array<double, 4> HGCalHistoSeedingImpl::boundaries() {
0556   switch (seedingSpace_) {
0557     case RPhi:
0558       return {{kROverZMin_, kROverZMax_, -M_PI, M_PI}};
0559     case XY:
0560       return {{-kXYMax_, kXYMax_, -kXYMax_, kXYMax_}};
0561   }
0562   return {{0., 0., 0., 0.}};
0563 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!