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)));
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);
0193
0194 content += histoClusters.at(z_side, bin1, binToSumRight).values[Bin::Content::Sum] /
0195 pow(2, bin22);
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.;
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
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
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
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
0521 Histogram histoCluster = fillHistoClusters(clustersPtrs);
0522
0523 Histogram smoothHistoCluster;
0524 if (seedingSpace_ == RPhi) {
0525
0526 Histogram smoothPhiHistoCluster = fillSmoothPhiHistoClusters(histoCluster, binsSumsHisto_);
0527
0528
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
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
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 }