File indexing completed on 2024-04-06 12:20:38
0001 #ifndef __L1Trigger_L1THGCal_HGCalHistoSeedingImpl_h__
0002 #define __L1Trigger_L1THGCal_HGCalHistoSeedingImpl_h__
0003
0004 #include "DataFormats/L1THGCal/interface/HGCalCluster.h"
0005 #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008
0009 #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h"
0010 #include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h"
0011 #include "L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h"
0012 #include "L1Trigger/L1THGCal/interface/backend/HGCalTriggerClusterIdentificationBase.h"
0013
0014 class HGCalHistoSeedingImpl {
0015 private:
0016 struct Bin {
0017 enum Content { Sum, Ecal, Hcal };
0018 std::array<float, 3> values = {{0., 0., 0.}};
0019 float weighted_x = 0.;
0020 float weighted_y = 0.;
0021 };
0022 template <typename T>
0023 class HistogramT {
0024 public:
0025 using Data = std::vector<T>;
0026 using iterator = typename Data::iterator;
0027 using const_iterator = typename Data::const_iterator;
0028
0029 public:
0030 HistogramT() : bins1_(0), bins2_(0), bins_(0) {}
0031 HistogramT(unsigned bins1, unsigned bins2)
0032 : bins1_(bins1), bins2_(bins2), bins_(bins1 * bins2), histogram_(bins_ * kSides_) {}
0033
0034 T& at(int zside, unsigned x1, unsigned x2) { return histogram_[index(zside, x1, x2)]; }
0035
0036 const T& at(int zside, unsigned x1, unsigned x2) const { return histogram_[index(zside, x1, x2)]; }
0037
0038 iterator begin() { return histogram_.begin(); }
0039 const_iterator begin() const { return histogram_.begin(); }
0040 iterator end() { return histogram_.end(); }
0041 const_iterator end() const { return histogram_.end(); }
0042
0043 private:
0044 static constexpr unsigned kSides_ = 2;
0045 unsigned bins1_ = 0;
0046 unsigned bins2_ = 0;
0047 unsigned bins_ = 0;
0048 Data histogram_;
0049
0050 unsigned index(int zside, unsigned x1, unsigned x2) const {
0051 if (x1 >= bins1_ || x2 >= bins2_) {
0052 throw cms::Exception("OutOfBound") << "Trying to access bin (" << x1 << "," << x2
0053 << ") in seeding histogram of size (" << bins1_ << "," << bins2_ << ")";
0054 }
0055 return x2 + bins2_ * x1 + bins_ * (zside > 0 ? 1 : 0);
0056 }
0057 };
0058 using Histogram = HistogramT<Bin>;
0059
0060 class Navigator {
0061 public:
0062 enum class AxisType { Bounded, Circular };
0063
0064 Navigator() : axis_types_({{AxisType::Bounded, AxisType::Bounded}}), bins_({{0, 0}}), home_({{0, 0}}) {}
0065
0066 Navigator(int bins1, AxisType type_axis1, int bins2, AxisType type_axis2)
0067 : axis_types_({{type_axis1, type_axis2}}), bins_({{bins1, bins2}}), home_({{0, 0}}) {}
0068
0069 void setHome(int x1, int x2) {
0070 if (x1 < 0 || x2 < 0 || x1 >= std::get<0>(bins_) || x2 >= std::get<1>(bins_)) {
0071 throw cms::Exception("OutOfBound") << "Setting invalid navigator home position (" << x1 << "," << x2 << "\n)";
0072 }
0073 home_[0] = x1;
0074 home_[1] = x2;
0075 }
0076
0077 std::array<int, 2> move(int offset1, int offset2) { return {{offset<0>(offset1), offset<1>(offset2)}}; }
0078
0079 template <unsigned N>
0080 int offset(int shift) {
0081 int shifted = std::get<N>(home_);
0082 int max = std::get<N>(bins_);
0083 switch (std::get<N>(axis_types_)) {
0084 case AxisType::Bounded:
0085 shifted += shift;
0086 if (shifted < 0 || shifted >= max)
0087 shifted = -1;
0088 break;
0089 case AxisType::Circular:
0090 shifted += shift;
0091 while (shifted < 0)
0092 shifted += max;
0093 while (shifted >= max)
0094 shifted -= max;
0095 break;
0096 default:
0097 break;
0098 }
0099 return shifted;
0100 }
0101
0102 private:
0103 std::array<AxisType, 2> axis_types_;
0104 std::array<int, 2> bins_;
0105 std::array<int, 2> home_;
0106 };
0107
0108 public:
0109 HGCalHistoSeedingImpl(const edm::ParameterSet& conf);
0110
0111 void setGeometry(const HGCalTriggerGeometryBase* const geom) { triggerTools_.setGeometry(geom); }
0112
0113 float dR(const l1t::HGCalCluster& clu, const GlobalPoint& seed) const;
0114
0115 void findHistoSeeds(const std::vector<edm::Ptr<l1t::HGCalCluster>>& clustersPtr,
0116 std::vector<std::pair<GlobalPoint, double>>& seedPositionsEnergy);
0117
0118 private:
0119 enum SeedingType { HistoMaxC3d, HistoSecondaryMaxC3d, HistoThresholdC3d, HistoInterpolatedMaxC3d };
0120 enum SeedingSpace { RPhi, XY };
0121 enum SeedingPosition { BinCentre, TCWeighted };
0122
0123 Histogram fillHistoClusters(const std::vector<edm::Ptr<l1t::HGCalCluster>>& clustersPtrs);
0124
0125 Histogram fillSmoothHistoClusters(const Histogram&, const vector<double>&, Bin::Content);
0126 Histogram fillSmoothPhiHistoClusters(const Histogram& histoClusters, const vector<unsigned>& binSums);
0127 Histogram fillSmoothRPhiHistoClusters(const Histogram& histoClusters);
0128
0129 void setSeedEnergyAndPosition(std::vector<std::pair<GlobalPoint, double>>& seedPositionsEnergy,
0130 int z_side,
0131 unsigned bin_R,
0132 unsigned bin_phi,
0133 const Bin& histBin);
0134
0135 std::vector<std::pair<GlobalPoint, double>> computeMaxSeeds(const Histogram& histoClusters);
0136
0137 std::vector<std::pair<GlobalPoint, double>> computeSecondaryMaxSeeds(const Histogram& histoClusters);
0138
0139 std::vector<std::pair<GlobalPoint, double>> computeInterpolatedMaxSeeds(const Histogram& histoClusters);
0140
0141 std::vector<std::pair<GlobalPoint, double>> computeThresholdSeeds(const Histogram& histoClusters);
0142
0143 std::array<double, 4> boundaries();
0144
0145 std::string seedingAlgoType_;
0146 SeedingType seedingType_;
0147 SeedingSpace seedingSpace_;
0148 SeedingPosition seedingPosition_;
0149
0150 unsigned nBins1_ = 42;
0151 unsigned nBins2_ = 216;
0152 std::vector<unsigned> binsSumsHisto_;
0153 double histoThreshold_ = 20.;
0154 static constexpr double area_per_triggercell_ =
0155 4.91E-05;
0156 std::vector<double> neighbour_weights_;
0157 std::vector<double> smoothing_ecal_;
0158 std::vector<double> smoothing_hcal_;
0159 bool seeds_norm_by_area_;
0160
0161 float area_10pct_;
0162
0163 HGCalTriggerTools triggerTools_;
0164 Navigator navigator_;
0165
0166 static constexpr unsigned neighbour_weights_size_ = 9;
0167 const double kROverZMin_ = 0.076;
0168 const double kROverZMax_ = 0.58;
0169
0170 static constexpr double kXYMax_ = 0.6;
0171 };
0172
0173 #endif
