Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:55:14 

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