CMS 3D CMS Logo
HGCalHistoSeedingImpl.cc
Go to the documentation of this file.
1 #include "L1Trigger/L1THGCal/interface/backend/HGCalHistoSeedingImpl.h"
2 #include "L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h"
3 #include "DataFormats/Math/interface/deltaR.h"
4 #include "DataFormats/Math/interface/deltaPhi.h"
5 #include <numeric>
6 
7 HGCalHistoSeedingImpl::HGCalHistoSeedingImpl(const edm::ParameterSet& conf)
8  : seedingAlgoType_(conf.getParameter<std::string>("type_histoalgo")),
9  nBins1_(conf.getParameter<unsigned>("nBins_X1_histo_multicluster")),
10  nBins2_(conf.getParameter<unsigned>("nBins_X2_histo_multicluster")),
11  binsSumsHisto_(conf.getParameter<std::vector<unsigned>>("binSumsHisto")),
12  histoThreshold_(conf.getParameter<double>("threshold_histo_multicluster")),
13  neighbour_weights_(conf.getParameter<std::vector<double>>("neighbour_weights")),
14  smoothing_ecal_(conf.getParameter<std::vector<double>>("seed_smoothing_ecal")),
15  smoothing_hcal_(conf.getParameter<std::vector<double>>("seed_smoothing_hcal")),
16  kROverZMin_(conf.getParameter<double>("kROverZMin")),
17  kROverZMax_(conf.getParameter<double>("kROverZMax")) {
18  if (seedingAlgoType_ == "HistoMaxC3d") {
19  seedingType_ = HistoMaxC3d;
20  } else if (seedingAlgoType_ == "HistoSecondaryMaxC3d") {
21  seedingType_ = HistoSecondaryMaxC3d;
22  } else if (seedingAlgoType_ == "HistoThresholdC3d") {
23  seedingType_ = HistoThresholdC3d;
24  } else if (seedingAlgoType_ == "HistoInterpolatedMaxC3d") {
25  seedingType_ = HistoInterpolatedMaxC3d;
26  } else {
27  throw cms::Exception("HGCTriggerParameterError") << "Unknown Multiclustering type '" << seedingAlgoType_;
28  }
29 
30  if (conf.getParameter<std::string>("seed_position") == "BinCentre") {
31  seedingPosition_ = BinCentre;
32  } else if (conf.getParameter<std::string>("seed_position") == "TCWeighted") {
33  seedingPosition_ = TCWeighted;
34  } else {
35  throw cms::Exception("HGCTriggerParameterError")
36  << "Unknown Seed Position option '" << conf.getParameter<std::string>("seed_position");
37  }
38  if (conf.getParameter<std::string>("seeding_space") == "RPhi") {
39  seedingSpace_ = RPhi;
40  navigator_ = Navigator(nBins1_, Navigator::AxisType::Bounded, nBins2_, Navigator::AxisType::Circular);
41  } else if (conf.getParameter<std::string>("seeding_space") == "XY") {
42  seedingSpace_ = XY;
43  navigator_ = Navigator(nBins1_, Navigator::AxisType::Bounded, nBins2_, Navigator::AxisType::Bounded);
44  } else {
45  throw cms::Exception("HGCTriggerParameterError")
46  << "Unknown seeding space '" << conf.getParameter<std::string>("seeding_space");
47  }
48 
49  edm::LogInfo("HGCalMulticlusterParameters")
50  << "\nMulticluster number of X1-bins for the histo algorithm: " << nBins1_
51  << "\nMulticluster number of X2-bins for the histo algorithm: " << nBins2_
52  << "\nMulticluster MIPT threshold for histo threshold algorithm: " << histoThreshold_
53  << "\nMulticluster type of multiclustering algortihm: " << seedingAlgoType_;
54 
55  if (seedingAlgoType_.find("Histo") != std::string::npos && seedingSpace_ == RPhi &&
56  nBins1_ != binsSumsHisto_.size()) {
57  throw cms::Exception("Inconsistent bin size")
58  << "Inconsistent nBins_X1_histo_multicluster ( " << nBins1_ << " ) and binSumsHisto ( " << binsSumsHisto_.size()
59  << " ) size in HGCalMulticlustering\n";
60  }
61 
62  if (neighbour_weights_.size() != neighbour_weights_size_) {
63  throw cms::Exception("Inconsistent vector size")
64  << "Inconsistent size of neighbour weights vector in HGCalMulticlustering ( " << neighbour_weights_.size()
65  << " ). Should be " << neighbour_weights_size_ << "\n";
66  }
67 }
68 
69 HGCalHistoSeedingImpl::Histogram HGCalHistoSeedingImpl::fillHistoClusters(
70  const std::vector<edm::Ptr<l1t::HGCalCluster>>& clustersPtrs) {
71  Histogram histoClusters(nBins1_, nBins2_);
72  std::array<double, 4> bounds = boundaries();
73  double minx1 = std::get<0>(bounds);
74  double maxx1 = std::get<1>(bounds);
75  double minx2 = std::get<2>(bounds);
76  double maxx2 = std::get<3>(bounds);
77 
78  for (auto& clu : clustersPtrs) {
79  float x1 = 0., x2 = 0;
80  switch (seedingSpace_) {
81  case RPhi:
82  x1 = sqrt(pow(clu->centreProj().x(), 2) + pow(clu->centreProj().y(), 2));
83  x2 = reco::reduceRange(clu->phi());
84  break;
85  case XY:
86  x1 = clu->centreProj().x();
87  x2 = clu->centreProj().y();
88  break;
89  };
90  if (x1 < minx1 || x1 >= maxx1) {
91  throw cms::Exception("OutOfBound") << "TC X1 = " << x1 << " out of the seeding histogram bounds " << minx1
92  << " - " << maxx1;
93  }
94  if (x2 < minx2 || x2 >= maxx2) {
95  throw cms::Exception("OutOfBound") << "TC X2 = " << x2 << " out of the seeding histogram bounds " << minx2
96  << " - " << maxx2;
97  }
98  unsigned bin1 = unsigned((x1 - minx1) * nBins1_ / (maxx1 - minx1));
99  unsigned bin2 = unsigned((x2 - minx2) * nBins2_ / (maxx2 - minx2));
100 
101  auto& bin = histoClusters.at(triggerTools_.zside(clu->detId()), bin1, bin2);
102  bin.values[Bin::Content::Sum] += clu->mipPt();
103  if (triggerTools_.isEm(clu->detId())) {
104  bin.values[Bin::Content::Ecal] += clu->mipPt();
105  } else {
106  bin.values[Bin::Content::Hcal] += clu->mipPt();
107  }
108  bin.weighted_x += (clu->centreProj().x()) * clu->mipPt();
109  bin.weighted_y += (clu->centreProj().y()) * clu->mipPt();
110  }
111 
112  for (auto& bin : histoClusters) {
113  bin.weighted_x /= bin.values[Bin::Content::Sum];
114  bin.weighted_y /= bin.values[Bin::Content::Sum];
115  }
116 
117  return histoClusters;
118 }
119 
120 HGCalHistoSeedingImpl::Histogram HGCalHistoSeedingImpl::fillSmoothHistoClusters(const Histogram& histoClusters,
121  const vector<double>& kernel,
122  Bin::Content binContent) {
123  Histogram histoSmooth(histoClusters);
124 
125  unsigned kernel_size = std::sqrt(kernel.size());
126  if (kernel_size * kernel_size != kernel.size()) {
127  throw cms::Exception("HGCTriggerParameterError") << "Only square kernels can be used.";
128  }
129  if (kernel_size % 2 != 1) {
130  throw cms::Exception("HGCTriggerParameterError") << "The kernel size must be an odd value.";
131  }
132  int shift_max = (kernel_size - 1) / 2;
133  double normalization = std::accumulate(kernel.begin(), kernel.end(), 0.);
134  for (int z_side : {-1, 1}) {
135  for (unsigned x1 = 0; x1 < nBins1_; x1++) {
136  for (unsigned x2 = 0; x2 < nBins2_; x2++) {
137  const auto& bin_orig = histoClusters.at(z_side, x1, x2);
138  double smooth = 0.;
139  navigator_.setHome(x1, x2);
140  for (int x1_shift = -shift_max; x1_shift <= shift_max; x1_shift++) {
141  int index1 = x1_shift + shift_max;
142  for (int x2_shift = -shift_max; x2_shift <= shift_max; x2_shift++) {
143  auto shifted = navigator_.move(x1_shift, x2_shift);
144  int index2 = x2_shift + shift_max;
145  double kernel_value = kernel.at(index1 * kernel_size + index2);
146  bool out = shifted[0] == -1 || shifted[1] == -1;
147  double content = (out ? 0. : histoClusters.at(z_side, shifted[0], shifted[1]).values[binContent]);
148  smooth += content * kernel_value;
149  }
150  }
151  auto& bin = histoSmooth.at(z_side, x1, x2);
152  bin.values[binContent] = smooth / normalization;
153  bin.weighted_x = bin_orig.weighted_x;
154  bin.weighted_y = bin_orig.weighted_y;
155  }
156  }
157  }
158 
159  return histoSmooth;
160 }
161 
162 HGCalHistoSeedingImpl::Histogram HGCalHistoSeedingImpl::fillSmoothPhiHistoClusters(const Histogram& histoClusters,
163  const vector<unsigned>& binSums) {
164  Histogram histoSumPhiClusters(nBins1_, nBins2_);
165 
166  for (int z_side : {-1, 1}) {
167  for (unsigned bin1 = 0; bin1 < nBins1_; bin1++) {
168  int nBinsSide = (binSums[bin1] - 1) / 2;
169  float R1 = kROverZMin_ + bin1 * (kROverZMax_ - kROverZMin_) / nBins1_;
170  float R2 = R1 + ((kROverZMax_ - kROverZMin_) / nBins1_);
171  double area =
172  ((M_PI * (pow(R2, 2) - pow(R1, 2))) / nBins2_) *
173  (1 +
174  2.0 *
175  (1 -
176  pow(0.5,
177  nBinsSide))); // Takes into account different area of bins in different R-rings + sum of quadratic weights used
178 
179  for (unsigned bin2 = 0; bin2 < nBins2_; bin2++) {
180  const auto& bin_orig = histoClusters.at(z_side, bin1, bin2);
181  float content = bin_orig.values[Bin::Content::Sum];
182 
183  for (int bin22 = 1; bin22 <= nBinsSide; bin22++) {
184  int binToSumLeft = bin2 - bin22;
185  if (binToSumLeft < 0)
186  binToSumLeft += nBins2_;
187  unsigned binToSumRight = bin2 + bin22;
188  if (binToSumRight >= nBins2_)
189  binToSumRight -= nBins2_;
190 
191  content += histoClusters.at(z_side, bin1, binToSumLeft).values[Bin::Content::Sum] /
192  pow(2, bin22); // quadratic kernel
193 
194  content += histoClusters.at(z_side, bin1, binToSumRight).values[Bin::Content::Sum] /
195  pow(2, bin22); // quadratic kernel
196  }
197 
198  auto& bin = histoSumPhiClusters.at(z_side, bin1, bin2);
199  bin.values[Bin::Content::Sum] = (content / area) * area_per_triggercell_;
200  bin.weighted_x = bin_orig.weighted_x;
201  bin.weighted_y = bin_orig.weighted_y;
202  }
203  }
204  }
205 
206  return histoSumPhiClusters;
207 }
208 
209 HGCalHistoSeedingImpl::Histogram HGCalHistoSeedingImpl::fillSmoothRPhiHistoClusters(const Histogram& histoClusters) {
210  Histogram histoSumRPhiClusters(nBins1_, nBins2_);
211 
212  for (int z_side : {-1, 1}) {
213  for (unsigned bin1 = 0; bin1 < nBins1_; bin1++) {
214  float weight =
215  (bin1 == 0 || bin1 == nBins1_ - 1) ? 1.5 : 2.; //Take into account edges with only one side up or down
216 
217  for (unsigned bin2 = 0; bin2 < nBins2_; bin2++) {
218  const auto& bin_orig = histoClusters.at(z_side, bin1, bin2);
219  float content = bin_orig.values[Bin::Content::Sum];
220  float contentDown = bin1 > 0 ? histoClusters.at(z_side, bin1 - 1, bin2).values[Bin::Content::Sum] : 0;
221  float contentUp = bin1 < (nBins1_ - 1) ? histoClusters.at(z_side, bin1 + 1, bin2).values[Bin::Content::Sum] : 0;
222 
223  auto& bin = histoSumRPhiClusters.at(z_side, bin1, bin2);
224  bin.values[Bin::Content::Sum] = (content + 0.5 * contentDown + 0.5 * contentUp) / weight;
225  bin.weighted_x = bin_orig.weighted_x;
226  bin.weighted_y = bin_orig.weighted_y;
227  }
228  }
229  }
230 
231  return histoSumRPhiClusters;
232 }
233 
234 void HGCalHistoSeedingImpl::setSeedEnergyAndPosition(std::vector<std::pair<GlobalPoint, double>>& seedPositionsEnergy,
235  int z_side,
236  unsigned bin1,
237  unsigned bin2,
238  const Bin& histBin) {
239  float x_seed = 0;
240  float y_seed = 0;
241  std::array<double, 4> bounds = boundaries();
242  double minx1 = std::get<0>(bounds);
243  double maxx1 = std::get<1>(bounds);
244  double minx2 = std::get<2>(bounds);
245  double maxx2 = std::get<3>(bounds);
246 
247  if (seedingPosition_ == BinCentre) {
248  float x1_seed = minx1 + (bin1 + 0.5) * (maxx1 - minx1) / nBins1_;
249  float x2_seed = minx2 + (bin2 + 0.5) * (maxx2 - minx2) / nBins2_;
250  switch (seedingSpace_) {
251  case RPhi:
252  x_seed = x1_seed * cos(x2_seed);
253  y_seed = x1_seed * sin(x2_seed);
254  break;
255  case XY:
256  x_seed = x1_seed;
257  y_seed = x2_seed;
258  };
259  } else if (seedingPosition_ == TCWeighted) {
260  x_seed = histBin.weighted_x;
261  y_seed = histBin.weighted_y;
262  }
263 
264  seedPositionsEnergy.emplace_back(GlobalPoint(x_seed, y_seed, z_side), histBin.values[Bin::Content::Sum]);
265 }
266 
267 std::vector<std::pair<GlobalPoint, double>> HGCalHistoSeedingImpl::computeMaxSeeds(const Histogram& histoClusters) {
268  std::vector<std::pair<GlobalPoint, double>> seedPositionsEnergy;
269 
270  for (int z_side : {-1, 1}) {
271  for (unsigned bin1 = 0; bin1 < nBins1_; bin1++) {
272  for (unsigned bin2 = 0; bin2 < nBins2_; bin2++) {
273  float MIPT_seed = histoClusters.at(z_side, bin1, bin2).values[Bin::Content::Sum];
274  bool isMax = MIPT_seed > histoThreshold_;
275  if (!isMax)
276  continue;
277 
278  navigator_.setHome(bin1, bin2);
279  auto pos_N = navigator_.move(1, 0);
280  auto pos_S = navigator_.move(-1, 0);
281  auto pos_W = navigator_.move(0, -1);
282  auto pos_E = navigator_.move(0, 1);
283  auto pos_NW = navigator_.move(1, -1);
284  auto pos_NE = navigator_.move(1, 1);
285  auto pos_SW = navigator_.move(-1, -1);
286  auto pos_SE = navigator_.move(-1, 1);
287 
288  float MIPT_N = (pos_N[0] != -1 && pos_N[1] != -1)
289  ? histoClusters.at(z_side, pos_N[0], pos_N[1]).values[Bin::Content::Sum]
290  : 0;
291  float MIPT_S = (pos_S[0] != -1 && pos_S[1] != -1)
292  ? histoClusters.at(z_side, pos_S[0], pos_S[1]).values[Bin::Content::Sum]
293  : 0;
294  float MIPT_W = (pos_W[0] != -1 && pos_W[1] != -1)
295  ? histoClusters.at(z_side, pos_W[0], pos_W[1]).values[Bin::Content::Sum]
296  : 0;
297  float MIPT_E = (pos_E[0] != -1 && pos_E[1] != -1)
298  ? histoClusters.at(z_side, pos_E[0], pos_E[1]).values[Bin::Content::Sum]
299  : 0;
300  float MIPT_NW = (pos_NW[0] != -1 && pos_NW[1] != -1)
301  ? histoClusters.at(z_side, pos_NW[0], pos_NW[1]).values[Bin::Content::Sum]
302  : 0;
303  float MIPT_NE = (pos_NE[0] != -1 && pos_NE[1] != -1)
304  ? histoClusters.at(z_side, pos_NE[0], pos_NE[1]).values[Bin::Content::Sum]
305  : 0;
306  float MIPT_SW = (pos_SW[0] != -1 && pos_SW[1] != -1)
307  ? histoClusters.at(z_side, pos_SW[0], pos_SW[1]).values[Bin::Content::Sum]
308  : 0;
309  float MIPT_SE = (pos_SE[0] != -1 && pos_SE[1] != -1)
310  ? histoClusters.at(z_side, pos_SE[0], pos_SE[1]).values[Bin::Content::Sum]
311  : 0;
312 
313  isMax &= MIPT_seed >= MIPT_S && MIPT_seed > MIPT_N && MIPT_seed >= MIPT_E && MIPT_seed >= MIPT_SE &&
314  MIPT_seed >= MIPT_NE && MIPT_seed > MIPT_W && MIPT_seed > MIPT_SW && MIPT_seed > MIPT_NW;
315 
316  if (isMax) {
317  setSeedEnergyAndPosition(seedPositionsEnergy, z_side, bin1, bin2, histoClusters.at(z_side, bin1, bin2));
318  }
319  }
320  }
321  }
322 
323  return seedPositionsEnergy;
324 }
325 
326 std::vector<std::pair<GlobalPoint, double>> HGCalHistoSeedingImpl::computeInterpolatedMaxSeeds(
327  const Histogram& histoClusters) {
328  std::vector<std::pair<GlobalPoint, double>> seedPositionsEnergy;
329 
330  for (int z_side : {-1, 1}) {
331  for (unsigned bin1 = 0; bin1 < nBins1_; bin1++) {
332  for (unsigned bin2 = 0; bin2 < nBins2_; bin2++) {
333  float MIPT_seed = histoClusters.at(z_side, bin1, bin2).values[Bin::Content::Sum];
334 
335  navigator_.setHome(bin1, bin2);
336  auto pos_N = navigator_.move(1, 0);
337  auto pos_S = navigator_.move(-1, 0);
338  auto pos_W = navigator_.move(0, -1);
339  auto pos_E = navigator_.move(0, 1);
340  auto pos_NW = navigator_.move(1, -1);
341  auto pos_NE = navigator_.move(1, 1);
342  auto pos_SW = navigator_.move(-1, -1);
343  auto pos_SE = navigator_.move(-1, 1);
344 
345  float MIPT_N = (pos_N[0] != -1 && pos_N[1] != -1)
346  ? histoClusters.at(z_side, pos_N[0], pos_N[1]).values[Bin::Content::Sum]
347  : 0;
348  float MIPT_S = (pos_S[0] != -1 && pos_S[1] != -1)
349  ? histoClusters.at(z_side, pos_S[0], pos_S[1]).values[Bin::Content::Sum]
350  : 0;
351  float MIPT_W = (pos_W[0] != -1 && pos_W[1] != -1)
352  ? histoClusters.at(z_side, pos_W[0], pos_W[1]).values[Bin::Content::Sum]
353  : 0;
354  float MIPT_E = (pos_E[0] != -1 && pos_E[1] != -1)
355  ? histoClusters.at(z_side, pos_E[0], pos_E[1]).values[Bin::Content::Sum]
356  : 0;
357  float MIPT_NW = (pos_NW[0] != -1 && pos_NW[1] != -1)
358  ? histoClusters.at(z_side, pos_NW[0], pos_NW[1]).values[Bin::Content::Sum]
359  : 0;
360  float MIPT_NE = (pos_NE[0] != -1 && pos_NE[1] != -1)
361  ? histoClusters.at(z_side, pos_NE[0], pos_NE[1]).values[Bin::Content::Sum]
362  : 0;
363  float MIPT_SW = (pos_SW[0] != -1 && pos_SW[1] != -1)
364  ? histoClusters.at(z_side, pos_SW[0], pos_SW[1]).values[Bin::Content::Sum]
365  : 0;
366  float MIPT_SE = (pos_SE[0] != -1 && pos_SE[1] != -1)
367  ? histoClusters.at(z_side, pos_SE[0], pos_SE[1]).values[Bin::Content::Sum]
368  : 0;
369 
370  float MIPT_pred = neighbour_weights_.at(0) * MIPT_NW + neighbour_weights_.at(1) * MIPT_N +
371  neighbour_weights_.at(2) * MIPT_NE + neighbour_weights_.at(3) * MIPT_W +
372  neighbour_weights_.at(5) * MIPT_E + neighbour_weights_.at(6) * MIPT_SW +
373  neighbour_weights_.at(7) * MIPT_S + neighbour_weights_.at(8) * MIPT_SE;
374 
375  bool isMax = MIPT_seed >= (MIPT_pred + histoThreshold_);
376 
377  if (isMax) {
378  setSeedEnergyAndPosition(seedPositionsEnergy, z_side, bin1, bin2, histoClusters.at(z_side, bin1, bin2));
379  }
380  }
381  }
382  }
383 
384  return seedPositionsEnergy;
385 }
386 
387 std::vector<std::pair<GlobalPoint, double>> HGCalHistoSeedingImpl::computeThresholdSeeds(
388  const Histogram& histoClusters) {
389  std::vector<std::pair<GlobalPoint, double>> seedPositionsEnergy;
390 
391  for (int z_side : {-1, 1}) {
392  for (unsigned bin1 = 0; bin1 < nBins1_; bin1++) {
393  for (unsigned bin2 = 0; bin2 < nBins2_; bin2++) {
394  float MIPT_seed = histoClusters.at(z_side, bin1, bin2).values[Bin::Content::Sum];
395  bool isSeed = MIPT_seed > histoThreshold_;
396 
397  if (isSeed) {
398  setSeedEnergyAndPosition(seedPositionsEnergy, z_side, bin1, bin2, histoClusters.at(z_side, bin1, bin2));
399  }
400  }
401  }
402  }
403 
404  return seedPositionsEnergy;
405 }
406 
407 std::vector<std::pair<GlobalPoint, double>> HGCalHistoSeedingImpl::computeSecondaryMaxSeeds(
408  const Histogram& histoClusters) {
409  std::vector<std::pair<GlobalPoint, double>> seedPositionsEnergy;
410 
411  HistogramT<uint8_t> primarySeedPositions(nBins1_, nBins2_);
412  HistogramT<uint8_t> vetoPositions(nBins1_, nBins2_);
413 
414  //Search for primary seeds
415  for (int z_side : {-1, 1}) {
416  for (unsigned bin1 = 0; bin1 < nBins1_; bin1++) {
417  for (unsigned bin2 = 0; bin2 < nBins2_; bin2++) {
418  float MIPT_seed = histoClusters.at(z_side, bin1, bin2).values[Bin::Content::Sum];
419  bool isMax = MIPT_seed > histoThreshold_;
420 
421  if (!isMax)
422  continue;
423 
424  float MIPT_S = bin1 < (nBins1_ - 1) ? histoClusters.at(z_side, bin1 + 1, bin2).values[Bin::Content::Sum] : 0;
425  float MIPT_N = bin1 > 0 ? histoClusters.at(z_side, bin1 - 1, bin2).values[Bin::Content::Sum] : 0;
426 
427  int binLeft = bin2 - 1;
428  if (binLeft < 0)
429  binLeft += nBins2_;
430  unsigned binRight = bin2 + 1;
431  if (binRight >= nBins2_)
432  binRight -= nBins2_;
433 
434  float MIPT_W = histoClusters.at(z_side, bin1, binLeft).values[Bin::Content::Sum];
435  float MIPT_E = histoClusters.at(z_side, bin1, binRight).values[Bin::Content::Sum];
436  float MIPT_NW = bin1 > 0 ? histoClusters.at(z_side, bin1 - 1, binLeft).values[Bin::Content::Sum] : 0;
437  float MIPT_NE = bin1 > 0 ? histoClusters.at(z_side, bin1 - 1, binRight).values[Bin::Content::Sum] : 0;
438  float MIPT_SW =
439  bin1 < (nBins1_ - 1) ? histoClusters.at(z_side, bin1 + 1, binLeft).values[Bin::Content::Sum] : 0;
440  float MIPT_SE =
441  bin1 < (nBins1_ - 1) ? histoClusters.at(z_side, bin1 + 1, binRight).values[Bin::Content::Sum] : 0;
442 
443  isMax &= MIPT_seed >= MIPT_S && MIPT_seed > MIPT_N && MIPT_seed >= MIPT_E && MIPT_seed >= MIPT_SE &&
444  MIPT_seed >= MIPT_NE && MIPT_seed > MIPT_W && MIPT_seed > MIPT_SW && MIPT_seed > MIPT_NW;
445 
446  if (isMax) {
447  setSeedEnergyAndPosition(seedPositionsEnergy, z_side, bin1, bin2, histoClusters.at(z_side, bin1, bin2));
448 
449  primarySeedPositions.at(z_side, bin1, bin2) = true;
450 
451  vetoPositions.at(z_side, bin1, binLeft) = true;
452  vetoPositions.at(z_side, bin1, binRight) = true;
453  if (bin1 > 0) {
454  vetoPositions.at(z_side, bin1 - 1, bin2) = true;
455  vetoPositions.at(z_side, bin1 - 1, binRight) = true;
456  vetoPositions.at(z_side, bin1 - 1, binLeft) = true;
457  }
458  if (bin1 < (nBins1_ - 1)) {
459  vetoPositions.at(z_side, bin1 + 1, bin2) = true;
460  vetoPositions.at(z_side, bin1 + 1, binRight) = true;
461  vetoPositions.at(z_side, bin1 + 1, binLeft) = true;
462  }
463  }
464  }
465  }
466  }
467 
468  //Search for secondary seeds
469 
470  for (int z_side : {-1, 1}) {
471  for (unsigned bin1 = 0; bin1 < nBins1_; bin1++) {
472  for (unsigned bin2 = 0; bin2 < nBins2_; bin2++) {
473  //Cannot be a secondary seed if already a primary seed, or adjacent to primary seed
474  if (primarySeedPositions.at(z_side, bin1, bin2) || vetoPositions.at(z_side, bin1, bin2))
475  continue;
476 
477  float MIPT_seed = histoClusters.at(z_side, bin1, bin2).values[Bin::Content::Sum];
478  bool isMax = MIPT_seed > histoThreshold_;
479 
480  float MIPT_S = bin1 < (nBins1_ - 1) ? histoClusters.at(z_side, bin1 + 1, bin2).values[Bin::Content::Sum] : 0;
481  float MIPT_N = bin1 > 0 ? histoClusters.at(z_side, bin1 - 1, bin2).values[Bin::Content::Sum] : 0;
482 
483  int binLeft = bin2 - 1;
484  if (binLeft < 0)
485  binLeft += nBins2_;
486  unsigned binRight = bin2 + 1;
487  if (binRight >= nBins2_)
488  binRight -= nBins2_;
489 
490  float MIPT_W = histoClusters.at(z_side, bin1, binLeft).values[Bin::Content::Sum];
491  float MIPT_E = histoClusters.at(z_side, bin1, binRight).values[Bin::Content::Sum];
492  float MIPT_NW = bin1 > 0 ? histoClusters.at(z_side, bin1 - 1, binLeft).values[Bin::Content::Sum] : 0;
493  float MIPT_NE = bin1 > 0 ? histoClusters.at(z_side, bin1 - 1, binRight).values[Bin::Content::Sum] : 0;
494  float MIPT_SW =
495  bin1 < (nBins1_ - 1) ? histoClusters.at(z_side, bin1 + 1, binLeft).values[Bin::Content::Sum] : 0;
496  float MIPT_SE =
497  bin1 < (nBins1_ - 1) ? histoClusters.at(z_side, bin1 + 1, binRight).values[Bin::Content::Sum] : 0;
498 
499  isMax &= (((bin1 < nBins1_ - 1) && vetoPositions.at(z_side, bin1 + 1, bin2)) or MIPT_seed >= MIPT_S) &&
500  (((bin1 > 0) && vetoPositions.at(z_side, bin1 - 1, bin2)) or MIPT_seed > MIPT_N) &&
501  ((vetoPositions.at(z_side, bin1, binRight)) or MIPT_seed >= MIPT_E) &&
502  (((bin1 < nBins1_ - 1) && vetoPositions.at(z_side, bin1 + 1, binRight)) or MIPT_seed >= MIPT_SE) &&
503  (((bin1 > 0) && vetoPositions.at(z_side, bin1 - 1, binRight)) or MIPT_seed >= MIPT_NE) &&
504  ((vetoPositions.at(z_side, bin1, binLeft)) or MIPT_seed > MIPT_W) &&
505  (((bin1 < nBins1_ - 1) && vetoPositions.at(z_side, bin1 + 1, binLeft)) or MIPT_seed > MIPT_SW) &&
506  (((bin1 > 0) && vetoPositions.at(z_side, bin1 - 1, binLeft)) or MIPT_seed > MIPT_NW);
507 
508  if (isMax) {
509  setSeedEnergyAndPosition(seedPositionsEnergy, z_side, bin1, bin2, histoClusters.at(z_side, bin1, bin2));
510  }
511  }
512  }
513  }
514 
515  return seedPositionsEnergy;
516 }
517 
518 void HGCalHistoSeedingImpl::findHistoSeeds(const std::vector<edm::Ptr<l1t::HGCalCluster>>& clustersPtrs,
519  std::vector<std::pair<GlobalPoint, double>>& seedPositionsEnergy) {
520  /* put clusters into an r/z x phi histogram */
521  Histogram histoCluster = fillHistoClusters(clustersPtrs);
522 
523  Histogram smoothHistoCluster;
524  if (seedingSpace_ == RPhi) {
525  /* smoothen along the phi direction + normalize each bin to same area */
526  Histogram smoothPhiHistoCluster = fillSmoothPhiHistoClusters(histoCluster, binsSumsHisto_);
527 
528  /* smoothen along the r/z direction */
529  smoothHistoCluster = fillSmoothRPhiHistoClusters(smoothPhiHistoCluster);
530  } else if (seedingSpace_ == XY) {
531  smoothHistoCluster = fillSmoothHistoClusters(histoCluster, smoothing_ecal_, Bin::Content::Ecal);
532  smoothHistoCluster = fillSmoothHistoClusters(smoothHistoCluster, smoothing_hcal_, Bin::Content::Hcal);
533  // Update sum with smoothen ECAL + HCAL
534  for (int z_side : {-1, 1}) {
535  for (unsigned x1 = 0; x1 < nBins1_; x1++) {
536  for (unsigned x2 = 0; x2 < nBins2_; x2++) {
537  auto& bin = smoothHistoCluster.at(z_side, x1, x2);
538  bin.values[Bin::Content::Sum] = bin.values[Bin::Content::Ecal] + bin.values[Bin::Content::Hcal];
539  }
540  }
541  }
542  }
543 
544  /* seeds determined with local maximum criteria */
545  if (seedingType_ == HistoMaxC3d)
546  seedPositionsEnergy = computeMaxSeeds(smoothHistoCluster);
547  else if (seedingType_ == HistoThresholdC3d)
548  seedPositionsEnergy = computeThresholdSeeds(smoothHistoCluster);
549  else if (seedingType_ == HistoInterpolatedMaxC3d)
550  seedPositionsEnergy = computeInterpolatedMaxSeeds(smoothHistoCluster);
551  else if (seedingType_ == HistoSecondaryMaxC3d)
552  seedPositionsEnergy = computeSecondaryMaxSeeds(smoothHistoCluster);
553 }
554 
555 std::array<double, 4> HGCalHistoSeedingImpl::boundaries() {
556  switch (seedingSpace_) {
557  case RPhi:
558  return {{kROverZMin_, kROverZMax_, -M_PI, M_PI}};
559  case XY:
560  return {{-kXYMax_, kXYMax_, -kXYMax_, kXYMax_}};
561  }
562  return {{0., 0., 0., 0.}};
563 }
PostProcessor_cff.normalization
normalization
Definition: PostProcessor_cff.py:9
HGCalHistoSeedingImpl::navigator_
Navigator navigator_
Definition: HGCalHistoSeedingImpl.h:161
HGCalHistoSeedingImpl::computeSecondaryMaxSeeds
std::vector< std::pair< GlobalPoint, double > > computeSecondaryMaxSeeds(const Histogram &histoClusters)
Definition: HGCalHistoSeedingImpl.cc:407
HGCalHistoSeedingImpl::kROverZMax_
const double kROverZMax_
Definition: HGCalHistoSeedingImpl.h:165
HGCalHistoSeedingImpl::Bin::weighted_y
float weighted_y
Definition: HGCalHistoSeedingImpl.h:20
HGCalHistoSeedingImpl::area_per_triggercell_
static constexpr double area_per_triggercell_
Definition: HGCalHistoSeedingImpl.h:154
HGCalHistoSeedingImpl::XY
Definition: HGCalHistoSeedingImpl.h:120
HGCalHistoSeedingImpl::nBins2_
unsigned nBins2_
Definition: HGCalHistoSeedingImpl.h:151
hgcalBackEndLayer2Producer_cfi.binSums
binSums
Definition: hgcalBackEndLayer2Producer_cfi.py:12
HGCalHistoSeedingImpl::kROverZMin_
const double kROverZMin_
Definition: HGCalHistoSeedingImpl.h:164
deltaPhi.h
testProducerWithPsetDescEmpty_cfi.x2
x2
Definition: testProducerWithPsetDescEmpty_cfi.py:28
HGCalHistoSeedingImpl::HistogramT::at
T & at(int zside, unsigned x1, unsigned x2)
Definition: HGCalHistoSeedingImpl.h:34
HGCalHistoSeedingImpl::RPhi
Definition: HGCalHistoSeedingImpl.h:120
HGCalHistoSeedingImpl::boundaries
std::array< double, 4 > boundaries()
Definition: HGCalHistoSeedingImpl.cc:555
HGCalHistoSeedingImpl::neighbour_weights_
std::vector< double > neighbour_weights_
Definition: HGCalHistoSeedingImpl.h:156
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
HGCalHistoSeedingImpl.h
HGCalHistoSeedingImpl::triggerTools_
HGCalTriggerTools triggerTools_
Definition: HGCalHistoSeedingImpl.h:160
HGCalHistoSeedingImpl::Navigator::setHome
void setHome(int x1, int x2)
Definition: HGCalHistoSeedingImpl.h:69
HGCalHistoSeedingImpl::fillHistoClusters
Histogram fillHistoClusters(const std::vector< edm::Ptr< l1t::HGCalCluster >> &clustersPtrs)
Definition: HGCalHistoSeedingImpl.cc:69
HGCalShowerShape.h
HGCalHistoSeedingImpl::TCWeighted
Definition: HGCalHistoSeedingImpl.h:121
HGCalHistoSeedingImpl::Navigator::move
std::array< int, 2 > move(int offset1, int offset2)
Definition: HGCalHistoSeedingImpl.h:77
HGCalHistoSeedingImpl::computeInterpolatedMaxSeeds
std::vector< std::pair< GlobalPoint, double > > computeInterpolatedMaxSeeds(const Histogram &histoClusters)
Definition: HGCalHistoSeedingImpl.cc:326
HGCalHistoSeedingImpl::HistoThresholdC3d
Definition: HGCalHistoSeedingImpl.h:119
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
testProducerWithPsetDescEmpty_cfi.x1
x1
Definition: testProducerWithPsetDescEmpty_cfi.py:33
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
HGCalHistoSeedingImpl::HistoInterpolatedMaxC3d
Definition: HGCalHistoSeedingImpl.h:119
HGCalHistoSeedingImpl::fillSmoothRPhiHistoClusters
Histogram fillSmoothRPhiHistoClusters(const Histogram &histoClusters)
Definition: HGCalHistoSeedingImpl.cc:209
HGCalHistoSeedingImpl::HistogramT
Definition: HGCalHistoSeedingImpl.h:23
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
HGCalHistoSeedingImpl::Bin
Definition: HGCalHistoSeedingImpl.h:16
HGCalHistoSeedingImpl::seedingSpace_
SeedingSpace seedingSpace_
Definition: HGCalHistoSeedingImpl.h:147
HGCalHistoSeedingImpl::neighbour_weights_size_
static constexpr unsigned neighbour_weights_size_
Definition: HGCalHistoSeedingImpl.h:163
Sum
Definition: SiPixelActionExecutor.h:17
HGCalTriggerTools::zside
int zside(const DetId &) const
Definition: HGCalTriggerTools.cc:210
GlobalPoint
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
HGCalHistoSeedingImpl::Navigator
Definition: HGCalHistoSeedingImpl.h:60
HGCalHistoSeedingImpl::Bin::Content
Content
Definition: HGCalHistoSeedingImpl.h:17
HGCalTriggerTools::isEm
bool isEm(const DetId &) const
Definition: HGCalTriggerTools.cc:151
HGCalHistoSeedingImpl::Bin::values
std::array< float, 3 > values
Definition: HGCalHistoSeedingImpl.h:18
Skims_PA_cff.content
content
Definition: Skims_PA_cff.py:19
edm::ParameterSet
Definition: ParameterSet.h:47
HGCalHistoSeedingImpl::Navigator::AxisType::Bounded
deltaR.h
HGCalHistoSeedingImpl::fillSmoothHistoClusters
Histogram fillSmoothHistoClusters(const Histogram &, const vector< double > &, Bin::Content)
Definition: HGCalHistoSeedingImpl.cc:120
HGCalHistoSeedingImpl::computeThresholdSeeds
std::vector< std::pair< GlobalPoint, double > > computeThresholdSeeds(const Histogram &histoClusters)
Definition: HGCalHistoSeedingImpl.cc:387
HGCalHistoSeedingImpl::kXYMax_
static constexpr double kXYMax_
Definition: HGCalHistoSeedingImpl.h:167
HGCalHistoSeedingImpl::histoThreshold_
double histoThreshold_
Definition: HGCalHistoSeedingImpl.h:153
HGCalHistoSeedingImpl::nBins1_
unsigned nBins1_
Definition: HGCalHistoSeedingImpl.h:150
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:49
trackerHitRTTI::vector
Definition: trackerHitRTTI.h:21
HGCalHistoSeedingImpl::fillSmoothPhiHistoClusters
Histogram fillSmoothPhiHistoClusters(const Histogram &histoClusters, const vector< unsigned > &binSums)
Definition: HGCalHistoSeedingImpl.cc:162
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
HGCalHistoSeedingImpl::setSeedEnergyAndPosition
void setSeedEnergyAndPosition(std::vector< std::pair< GlobalPoint, double >> &seedPositionsEnergy, int z_side, unsigned bin_R, unsigned bin_phi, const Bin &histBin)
Definition: HGCalHistoSeedingImpl.cc:234
edm::Ptr
Definition: AssociationVector.h:31
HGCalHistoSeedingImpl::seedingAlgoType_
std::string seedingAlgoType_
Definition: HGCalHistoSeedingImpl.h:145
HGCalHistoSeedingImpl::seedingPosition_
SeedingPosition seedingPosition_
Definition: HGCalHistoSeedingImpl.h:148
newFWLiteAna.bin
bin
Definition: newFWLiteAna.py:161
std
Definition: JetResolutionObject.h:76
HGCalHistoSeedingImpl::HGCalHistoSeedingImpl
HGCalHistoSeedingImpl(const edm::ParameterSet &conf)
Definition: HGCalHistoSeedingImpl.cc:7
Exception
Definition: hltDiff.cc:245
HGCalHistoSeedingImpl::smoothing_hcal_
std::vector< double > smoothing_hcal_
Definition: HGCalHistoSeedingImpl.h:158
or
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
HGCalHistoSeedingImpl::binsSumsHisto_
std::vector< unsigned > binsSumsHisto_
Definition: HGCalHistoSeedingImpl.h:152
HGCalHistoSeedingImpl::computeMaxSeeds
std::vector< std::pair< GlobalPoint, double > > computeMaxSeeds(const Histogram &histoClusters)
Definition: HGCalHistoSeedingImpl.cc:267
HGCalHistoSeedingImpl::Bin::weighted_x
float weighted_x
Definition: HGCalHistoSeedingImpl.h:19
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
HGCalHistoSeedingImpl::HistoMaxC3d
Definition: HGCalHistoSeedingImpl.h:119
HGCalHistoSeedingImpl::seedingType_
SeedingType seedingType_
Definition: HGCalHistoSeedingImpl.h:146
HGCalHistoSeedingImpl::Navigator::AxisType::Circular
HGCalHistoSeedingImpl::findHistoSeeds
void findHistoSeeds(const std::vector< edm::Ptr< l1t::HGCalCluster >> &clustersPtr, std::vector< std::pair< GlobalPoint, double >> &seedPositionsEnergy)
Definition: HGCalHistoSeedingImpl.cc:518
HGCalHistoSeedingImpl::smoothing_ecal_
std::vector< double > smoothing_ecal_
Definition: HGCalHistoSeedingImpl.h:157
custom_jme_cff.area
area
Definition: custom_jme_cff.py:140
weight
Definition: weight.py:1
HGCalHistoSeedingImpl::HistoSecondaryMaxC3d
Definition: HGCalHistoSeedingImpl.h:119
HGCalHistoSeedingImpl::BinCentre
Definition: HGCalHistoSeedingImpl.h:121
reco::reduceRange
constexpr T reduceRange(T x)
Definition: deltaPhi.h:18
Generated for CMSSW Reference Manual by doxygen 1.8.16