FTLRecHit/interface/FTLCluster.h

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(); };
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(); };
0127     return weighted_mean(this->theHitENERGY, y_pos);
0128   }
0129
0130   inline float positionError(const float sigmaPos) const {
0131     float sumW2(0.f), sumW(0.f);
0132     for (const auto& hitW : theHitENERGY) {
0133       sumW2 += hitW * hitW;
0134       sumW += hitW;
0135     }
0136     if (sumW > 0)
0137       return sigmaPos * std::sqrt(sumW2) / sumW;
0138     else
0139       return -999.f;
0140   }
0141
0142   inline float time() const {
0143     auto t = [this](unsigned int i) { return this->theHitTIME[i]; };
0144     return weighted_mean(this->theHitENERGY, t);
0145   }
0146
0147   inline float timeError() const {
0148     auto t_err = [this](unsigned int i) { return this->theHitTIME_ERROR[i]; };
0149     return weighted_mean_error(this->theHitENERGY, t_err);
0150   }
0151
0152   // Return number of hits.
0153   inline int size() const { return theHitENERGY.size(); }
0154
0155   // Return cluster dimension in the x direction.
0156   inline int sizeX() const { return rowSpan() + 1; }
0157
0158   // Return cluster dimension in the y direction.
0159   inline int sizeY() const { return colSpan() + 1; }
0160
0161   inline float energy() const {
0162     return std::accumulate(theHitENERGY.begin(), theHitENERGY.end(), 0.f);
0163   }  // Return total cluster energy.
0164
0165   inline int minHitRow() const { return theMinHitRow; }             // The min x index.
0166   inline int maxHitRow() const { return minHitRow() + rowSpan(); }  // The max x index.
0167   inline int minHitCol() const { return theMinHitCol; }             // The min y index.
0168   inline int maxHitCol() const { return minHitCol() + colSpan(); }  // The max y index.
0169
0170   const std::vector<uint8_t>& hitOffset() const { return theHitOffset; }
0171   const std::vector<float>& hitENERGY() const { return theHitENERGY; }
0172   const std::vector<float>& hitTIME() const { return theHitTIME; }
0173   const std::vector<float>& hitTIME_ERROR() const { return theHitTIME_ERROR; }
0174
0175   // infinite faster than above...
0176   FTLHit hit(int i) const {
0177     return FTLHit(minHitRow() + theHitOffset[i * 2],
0178                   minHitCol() + theHitOffset[i * 2 + 1],
0179                   theHitENERGY[i],
0180                   theHitTIME[i],
0181                   theHitTIME_ERROR[i]);
0182   }
0183
0184   FTLHit seed() const { return hit(seed_); }
0185
0186   int colSpan() const { return theHitColSpan; }
0187
0188   int rowSpan() const { return theHitRowSpan; }
0189
0190   const DetId& id() const { return theid; }
0191   const DetId& detid() const { return id(); }
0192
0193   bool overflowCol() const { return overflow_(theHitColSpan); }
0194
0195   bool overflowRow() const { return overflow_(theHitRowSpan); }
0196
0197   bool overflow() const { return overflowCol() || overflowRow(); }
0198
0199   void packCol(uint16_t ymin, uint16_t yspan) {
0200     theMinHitCol = ymin;
0201     theHitColSpan = std::min(yspan, uint16_t(MAXSPAN));
0202   }
0203   void packRow(uint16_t xmin, uint16_t xspan) {
0204     theMinHitRow = xmin;
0205     theHitRowSpan = std::min(xspan, uint16_t(MAXSPAN));
0206   }
0207
0208   void setClusterPosX(float posx) { pos_x = posx; }
0209   void setClusterErrorX(float errx) { err_x = errx; }
0210   void setClusterErrorTime(float errtime) { err_time = errtime; }
0211   float getClusterPosX() const { return pos_x; }
0212   float getClusterErrorX() const { return err_x; }
0213   float getClusterErrorTime() const { return err_time; }
0214
0215 private:
0216   DetId theid;
0217
0218   std::vector<uint8_t> theHitOffset;
0219   std::vector<float> theHitENERGY;
0220   std::vector<float> theHitTIME;
0221   std::vector<float> theHitTIME_ERROR;
0222
0223   uint16_t theMinHitRow = MAXPOS;  // Minimum hit index in the x direction (low edge).
0224   uint16_t theMinHitCol = MAXPOS;  // Minimum hit index in the y direction (left edge).
0225   uint8_t theHitRowSpan = 0;       // Span hit index in the x direction (low edge).
0226   uint8_t theHitColSpan = 0;       // Span hit index in the y direction (left edge).
0227
0228   float pos_x = -99999.9f;  // For pixels with internal position information in one coordinate (i.e. BTL crystals)
0229   float err_x = -99999.9f;  // For pixels with internal position information in one coordinate (i.e. BTL crystals)
0230   float err_time = -99999.9f;
0231
0232   uint8_t seed_;
0233
0234   template <typename SumFunc, typename OutFunc>
0235   float weighted_sum(const std::vector<float>& weights, SumFunc&& sumFunc, OutFunc&& outFunc) const {
0236     float tot = 0;
0237     float sumW = 0;
0238     for (unsigned int i = 0; i < weights.size(); ++i) {
0239       tot += sumFunc(i);
0240       sumW += weights[i];
0241     }
0242     return outFunc(tot, sumW);
0243   }
0244
0245   template <typename Value>
0246   float weighted_mean(const std::vector<float>& weights, Value&& value) const {
0247     auto sumFunc = [&weights, value](unsigned int i) { return weights[i] * value(i); };
0248     auto outFunc = [](float x, float y) {
0249       if (y > 0)
0250         return (float)x / y;
0251       else
0252         return -999.f;
0253     };
0254     return weighted_sum(weights, sumFunc, outFunc);
0255   }
0256
0257   template <typename Err>
0258   float weighted_mean_error(const std::vector<float>& weights, Err&& err) const {
0259     auto sumFunc = [&weights, err](unsigned int i) { return weights[i] * weights[i] * err(i) * err(i); };
0260     auto outFunc = [](float x, float y) {
0261       if (y > 0)
0262         return (float)sqrt(x) / y;
0263       else
0264         return -999.f;
0265     };
0266     return weighted_sum(weights, sumFunc, outFunc);
0267   }
0268
0269   static int overflow_(uint16_t span) { return span == uint16_t(MAXSPAN); }
0270 };
0271
0272 // Comparison operators  (needed by DetSetVector & SortedCollection )
0273 inline bool operator<(const FTLCluster& one, const FTLCluster& other) {
0274   if (one.detid() == other.detid()) {
0275     if (one.minHitRow() < other.minHitRow()) {
0276       return true;
0277     } else if (one.minHitRow() > other.minHitRow()) {
0278       return false;
0279     } else if (one.minHitCol() < other.minHitCol()) {
0280       return true;
0281     } else {
0282       return false;
0283     }
0284   }
0285   return one.detid() < other.detid();
0286 }
0287
0288 inline bool operator<(const FTLCluster& one, const uint32_t& detid) { return one.detid() < detid; }
0289
0290 inline bool operator<(const uint32_t& detid, const FTLCluster& other) { return detid < other.detid(); }
0291
0292 #endif