Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-08 22:26:30

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;  // Hex_Wafer_Area (x/z units)/N_TC (per wafer) = (0.866*((hexWafer_minimal_diameter)*(1./319.))^2 / 48)
0156   std::vector<double> neighbour_weights_;
0157   std::vector<double> smoothing_ecal_;
0158   std::vector<double> smoothing_hcal_;
0159 
0160   HGCalTriggerTools triggerTools_;
0161   Navigator navigator_;
0162 
0163   static constexpr unsigned neighbour_weights_size_ = 9;
0164   const double kROverZMin_ = 0.076;
0165   const double kROverZMax_ = 0.58;
0166 
0167   static constexpr double kXYMax_ = 0.6;
0168 };
0169 
0170 #endif