Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-18 09:14:52

0001 #ifndef DataFormats_FTLRecHit_FTLCluster_h
0002 #define DataFormats_FTLRecHit_FTLCluster_h
0003 
0004 /** \class FTLCluster
0005  *  
0006  * based on SiPixelCluster
0007  *
0008  * \author Paolo Meridiani
0009  */
0010 
0011 #include <cmath>
0012 #include <vector>
0013 #include <cstdint>
0014 #include <cassert>
0015 #include <algorithm>
0016 #include <numeric>
0017 
0018 #include "DataFormats/DetId/interface/DetId.h"
0019 
0020 class FTLCluster {
0021 public:
0022   typedef DetId key_type;
0023 
0024   class FTLHit {
0025   public:
0026     constexpr FTLHit() : x_(0), y_(0), energy_(0), time_(0), time_error_(0) {}
0027     constexpr FTLHit(uint16_t hit_x, uint16_t hit_y, float hit_energy, float hit_time, float hit_time_error)
0028         : x_(hit_x), y_(hit_y), energy_(hit_energy), time_(hit_time), time_error_(hit_time_error) {}
0029     constexpr uint16_t x() { return x_; }
0030     constexpr uint16_t y() { return y_; }
0031     constexpr uint16_t energy() { return energy_; }
0032     constexpr uint16_t time() { return time_; }
0033     constexpr uint16_t time_error() { return time_error_; }
0034 
0035   private:
0036     uint16_t x_;  //row
0037     uint16_t y_;  //col
0038     float energy_;
0039     float time_;
0040     float time_error_;
0041   };
0042 
0043   //--- Integer shift in x and y directions.
0044   class Shift {
0045   public:
0046     constexpr Shift(int dx, int dy) : dx_(dx), dy_(dy) {}
0047     constexpr Shift() : dx_(0), dy_(0) {}
0048     constexpr int dx() const { return dx_; }
0049     constexpr int dy() const { return dy_; }
0050 
0051   private:
0052     int dx_;
0053     int dy_;
0054   };
0055 
0056   //--- Position of a FTL Hit
0057   class FTLHitPos {
0058   public:
0059     constexpr FTLHitPos() : row_(0), col_(0) {}
0060     constexpr FTLHitPos(int row, int col) : row_(row), col_(col) {}
0061     constexpr int row() const { return row_; }
0062     constexpr int col() const { return col_; }
0063     constexpr FTLHitPos operator+(const Shift& shift) const {
0064       return FTLHitPos(row() + shift.dx(), col() + shift.dy());
0065     }
0066 
0067   private:
0068     int row_;
0069     int col_;
0070   };
0071 
0072   static constexpr unsigned int MAXSPAN = 255;
0073   static constexpr unsigned int MAXPOS = 2047;
0074 
0075   /** Construct from a range of digis that form a cluster and from 
0076    *  a DetID. The range is assumed to be non-empty.
0077    */
0078   FTLCluster() {}
0079 
0080   FTLCluster(DetId id,
0081              unsigned int isize,
0082              float const* energys,
0083              float const* times,
0084              float const* time_errors,
0085              uint16_t const* xpos,
0086              uint16_t const* ypos,
0087              uint16_t const xmin,
0088              uint16_t const ymin)
0089       : theid(id),
0090         theHitOffset(2 * isize),
0091         theHitENERGY(energys, energys + isize),
0092         theHitTIME(times, times + isize),
0093         theHitTIME_ERROR(time_errors, time_errors + isize) {
0094     uint16_t maxCol = 0;
0095     uint16_t maxRow = 0;
0096     int maxHit = -1;
0097     float maxEnergy = -99999;
0098     for (unsigned int i = 0; i != isize; ++i) {
0099       uint16_t xoffset = xpos[i] - xmin;
0100       uint16_t yoffset = ypos[i] - ymin;
0101       theHitOffset[i * 2] = std::min(uint16_t(MAXSPAN), xoffset);
0102       theHitOffset[i * 2 + 1] = std::min(uint16_t(MAXSPAN), yoffset);
0103       if (xoffset > maxRow)
0104         maxRow = xoffset;
0105       if (yoffset > maxCol)
0106         maxCol = yoffset;
0107       if (theHitENERGY[i] > maxEnergy) {
0108         maxHit = i;
0109         maxEnergy = theHitENERGY[i];
0110       }
0111     }
0112     packRow(xmin, maxRow);
0113     packCol(ymin, maxCol);
0114 
0115     if (maxHit >= 0)
0116       seed_ = std::min(uint8_t(MAXSPAN), uint8_t(maxHit));
0117   }
0118 
0119   // linear average position (barycenter)
0120   inline float x() const {
0121     auto x_pos = [this](unsigned int i) { return this->theHitOffset[i * 2] + minHitRow() + 0.5f; };
0122     return weighted_mean(this->theHitENERGY, x_pos);
0123   }
0124 
0125   inline float y() const {
0126     auto y_pos = [this](unsigned int i) { return this->theHitOffset[i * 2 + 1] + minHitCol() + 0.5f; };
0127     return weighted_mean(this->theHitENERGY, y_pos);
0128   }
0129 
0130   inline float time() const {
0131     auto t = [this](unsigned int i) { return this->theHitTIME[i]; };
0132     return weighted_mean(this->theHitENERGY, t);
0133   }
0134 
0135   inline float timeError() const {
0136     auto t_err = [this](unsigned int i) { return this->theHitTIME_ERROR[i]; };
0137     return weighted_mean_error(this->theHitENERGY, t_err);
0138   }
0139 
0140   // Return number of hits.
0141   inline int size() const { return theHitENERGY.size(); }
0142 
0143   // Return cluster dimension in the x direction.
0144   inline int sizeX() const { return rowSpan() + 1; }
0145 
0146   // Return cluster dimension in the y direction.
0147   inline int sizeY() const { return colSpan() + 1; }
0148 
0149   inline float energy() const {
0150     return std::accumulate(theHitENERGY.begin(), theHitENERGY.end(), 0.f);
0151   }  // Return total cluster energy.
0152 
0153   inline int minHitRow() const { return theMinHitRow; }             // The min x index.
0154   inline int maxHitRow() const { return minHitRow() + rowSpan(); }  // The max x index.
0155   inline int minHitCol() const { return theMinHitCol; }             // The min y index.
0156   inline int maxHitCol() const { return minHitCol() + colSpan(); }  // The max y index.
0157 
0158   const std::vector<uint8_t>& hitOffset() const { return theHitOffset; }
0159   const std::vector<float>& hitENERGY() const { return theHitENERGY; }
0160   const std::vector<float>& hitTIME() const { return theHitTIME; }
0161   const std::vector<float>& hitTIME_ERROR() const { return theHitTIME_ERROR; }
0162 
0163   // infinite faster than above...
0164   FTLHit hit(int i) const {
0165     return FTLHit(minHitRow() + theHitOffset[i * 2],
0166                   minHitCol() + theHitOffset[i * 2 + 1],
0167                   theHitENERGY[i],
0168                   theHitTIME[i],
0169                   theHitTIME_ERROR[i]);
0170   }
0171 
0172   FTLHit seed() const { return hit(seed_); }
0173 
0174   int colSpan() const { return theHitColSpan; }
0175 
0176   int rowSpan() const { return theHitRowSpan; }
0177 
0178   const DetId& id() const { return theid; }
0179   const DetId& detid() const { return id(); }
0180 
0181   bool overflowCol() const { return overflow_(theHitColSpan); }
0182 
0183   bool overflowRow() const { return overflow_(theHitRowSpan); }
0184 
0185   bool overflow() const { return overflowCol() || overflowRow(); }
0186 
0187   void packCol(uint16_t ymin, uint16_t yspan) {
0188     theMinHitCol = ymin;
0189     theHitColSpan = std::min(yspan, uint16_t(MAXSPAN));
0190   }
0191   void packRow(uint16_t xmin, uint16_t xspan) {
0192     theMinHitRow = xmin;
0193     theHitRowSpan = std::min(xspan, uint16_t(MAXSPAN));
0194   }
0195 
0196   void setClusterErrorX(float errx) { err_x = errx; }
0197   void setClusterErrorY(float erry) { err_y = erry; }
0198   void setClusterErrorTime(float errtime) { err_time = errtime; }
0199   float getClusterErrorX() const { return err_x; }
0200   float getClusterErrorY() const { return err_y; }
0201   float getClusterErrorTime() const { return err_time; }
0202 
0203 private:
0204   DetId theid;
0205 
0206   std::vector<uint8_t> theHitOffset;
0207   std::vector<float> theHitENERGY;
0208   std::vector<float> theHitTIME;
0209   std::vector<float> theHitTIME_ERROR;
0210 
0211   uint16_t theMinHitRow = MAXPOS;  // Minimum hit index in the x direction (low edge).
0212   uint16_t theMinHitCol = MAXPOS;  // Minimum hit index in the y direction (left edge).
0213   uint8_t theHitRowSpan = 0;       // Span hit index in the x direction (low edge).
0214   uint8_t theHitColSpan = 0;       // Span hit index in the y direction (left edge).
0215 
0216   float err_x = -99999.9f;
0217   float err_y = -99999.9f;
0218   float err_time = -99999.9f;
0219 
0220   uint8_t seed_;
0221 
0222   template <typename SumFunc, typename OutFunc>
0223   float weighted_sum(const std::vector<float>& weights, SumFunc&& sumFunc, OutFunc&& outFunc) const {
0224     float tot = 0;
0225     float sumW = 0;
0226     for (unsigned int i = 0; i < weights.size(); ++i) {
0227       tot += sumFunc(i);
0228       sumW += weights[i];
0229     }
0230     return outFunc(tot, sumW);
0231   }
0232 
0233   template <typename Value>
0234   float weighted_mean(const std::vector<float>& weights, Value&& value) const {
0235     auto sumFunc = [&weights, value](unsigned int i) { return weights[i] * value(i); };
0236     auto outFunc = [](float x, float y) {
0237       if (y > 0)
0238         return (float)x / y;
0239       else
0240         return -999.f;
0241     };
0242     return weighted_sum(weights, sumFunc, outFunc);
0243   }
0244 
0245   template <typename Err>
0246   float weighted_mean_error(const std::vector<float>& weights, Err&& err) const {
0247     auto sumFunc = [&weights, err](unsigned int i) { return weights[i] * weights[i] * err(i) * err(i); };
0248     auto outFunc = [](float x, float y) {
0249       if (y > 0)
0250         return (float)sqrt(x) / y;
0251       else
0252         return -999.f;
0253     };
0254     return weighted_sum(weights, sumFunc, outFunc);
0255   }
0256 
0257   static int overflow_(uint16_t span) { return span == uint16_t(MAXSPAN); }
0258 };
0259 
0260 // Comparison operators  (needed by DetSetVector & SortedCollection )
0261 inline bool operator<(const FTLCluster& one, const FTLCluster& other) {
0262   if (one.detid() == other.detid()) {
0263     if (one.minHitRow() < other.minHitRow()) {
0264       return true;
0265     } else if (one.minHitRow() > other.minHitRow()) {
0266       return false;
0267     } else if (one.minHitCol() < other.minHitCol()) {
0268       return true;
0269     } else {
0270       return false;
0271     }
0272   }
0273   return one.detid() < other.detid();
0274 }
0275 
0276 inline bool operator<(const FTLCluster& one, const uint32_t& detid) { return one.detid() < detid; }
0277 
0278 inline bool operator<(const uint32_t& detid, const FTLCluster& other) { return detid < other.detid(); }
0279 
0280 #endif