FTLCluster.h

Go to the documentation of this file.

#ifndef DataFormats_FTLRecHit_FTLCluster_h
#define DataFormats_FTLRecHit_FTLCluster_h

#include <cmath>
#include <vector>
#include <cstdint>
#include <cassert>
#include <algorithm>
#include <numeric>
#include <functional>

#include "DataFormats/DetId/interface/DetId.h"

class FTLCluster {
public:

  typedef DetId key_type;

  class FTLHit {
  public:
    constexpr FTLHit() : x_(0), y_(0), energy_(0),  time_(0), time_error_(0) {}
    constexpr FTLHit(uint16_t hit_x, uint16_t hit_y, float hit_energy, float hit_time, float hit_time_error) :
    x_(hit_x), y_(hit_y), energy_(hit_energy), time_(hit_time), time_error_(hit_time_error) {}
    constexpr uint16_t x() { return x_; }
    constexpr uint16_t y() { return y_; }
    constexpr uint16_t energy() { return energy_; }
    constexpr uint16_t time() { return time_; }
    constexpr uint16_t time_error() { return time_error_; }
  private:
    uint16_t x_; //row
    uint16_t y_; //col
    float energy_;
    float time_;
    float time_error_;
  };
  
  //--- Integer shift in x and y directions.
  class Shift {
  public:
    constexpr Shift( int dx, int dy) : dx_(dx), dy_(dy) {}
    constexpr Shift() : dx_(0), dy_(0) {}
    constexpr int dx() const { return dx_;}
    constexpr int dy() const { return dy_;}
  private:
    int dx_;
    int dy_;
  };
  
  //--- Position of a FTL Hit
  class FTLHitPos {
  public:
    constexpr FTLHitPos() : row_(0), col_(0) {}
    constexpr FTLHitPos(int row, int col) : row_(row) , col_(col) {}
    constexpr int row() const { return row_;}
    constexpr int col() const { return col_;}
    constexpr FTLHitPos operator+( const Shift& shift) const {
      return FTLHitPos( row() + shift.dx(), col() + shift.dy());
    }
  private:
    int row_;
    int col_;
  };
    
  static constexpr unsigned int MAXSPAN=255;
  static constexpr unsigned int MAXPOS=2047;
  
  FTLCluster() {}
  
  FTLCluster(DetId id, unsigned int isize, float const * energys, float const* times, float const* time_errors,
        uint16_t const * xpos,  uint16_t const * ypos, 
        uint16_t const  xmin,  uint16_t const  ymin) :   
  theid(id), theHitOffset(2*isize), theHitENERGY(energys,energys+isize), theHitTIME(times,times+isize), theHitTIME_ERROR(time_errors,time_errors+isize)  {
    uint16_t maxCol = 0;
    uint16_t maxRow = 0;
    int maxHit=-1;
    float maxEnergy=-99999;
    for (unsigned int i=0; i!=isize; ++i) {
      uint16_t xoffset = xpos[i]-xmin;
      uint16_t yoffset = ypos[i]-ymin;
      theHitOffset[i*2] = std::min(uint16_t(MAXSPAN),xoffset);
      theHitOffset[i*2+1] = std::min(uint16_t(MAXSPAN),yoffset);
      if (xoffset > maxRow) maxRow = xoffset; 
      if (yoffset > maxCol) maxCol = yoffset; 
      if (theHitENERGY[i]>maxEnergy)
    {
      maxHit=i;
      maxEnergy=theHitENERGY[i];
    }
    }
    packRow(xmin,maxRow);
    packCol(ymin,maxCol);

    if (maxHit>=0)
      seed_=std::min(uint8_t(MAXSPAN),uint8_t(maxHit));
  }
  
  // linear average position (barycenter) 
  inline float x() const {
    auto x_pos=[this](unsigned int i) { return this->theHitOffset[i*2] + minHitRow() + 0.5f; };
    return weighted_mean(this->theHitENERGY,x_pos);
  }
  
  inline float y() const {
    auto y_pos=[this](unsigned int i) { return this->theHitOffset[i*2+1] + minHitCol() + 0.5f; };
    return weighted_mean(this->theHitENERGY,y_pos);
  }

  inline float time() const {
    auto t=[this](unsigned int i) { return this->theHitTIME[i]; };
    return weighted_mean(this->theHitENERGY,t);
  }

  inline float timeError() const {
    auto t_err=[this](unsigned int i) { return this->theHitTIME_ERROR[i]; };
    return weighted_mean_error(this->theHitENERGY,t_err);
  }
  
  // Return number of hits.
  inline int size() const { return theHitENERGY.size();}
  
  // Return cluster dimension in the x direction.
  inline int sizeX() const { return rowSpan() +1;}
  
  // Return cluster dimension in the y direction.
  inline int sizeY() const { return colSpan() +1;}
  
  
  inline float energy() const {
    return std::accumulate(theHitENERGY.begin(), theHitENERGY.end(),0.f);
  } // Return total cluster energy.
  
  inline int minHitRow() const { return theMinHitRow;} // The min x index.
  inline int maxHitRow() const { return minHitRow() + rowSpan();} // The max x index.
  inline int minHitCol() const { return theMinHitCol;} // The min y index.
  inline int maxHitCol() const { return minHitCol() + colSpan();} // The max y index.
  
  const std::vector<uint8_t> & hitOffset() const { return theHitOffset;}
  const std::vector<float> & hitENERGY() const { return theHitENERGY;}
  const std::vector<float> & hitTIME() const { return theHitTIME;}
  const std::vector<float> & hitTIME_ERROR() const { return theHitTIME_ERROR;}
    
  // infinite faster than above...
  FTLHit hit(int i) const {
    return FTLHit(minHitRow() + theHitOffset[i*2],
         minHitCol() + theHitOffset[i*2+1],
          theHitENERGY[i],
          theHitTIME[i],
          theHitTIME_ERROR[i]
         );

  }
  
  FTLHit seed() const {
    return hit(seed_);
  }
  
  int colSpan() const {return theHitColSpan; }
  
  int rowSpan() const { return theHitRowSpan; }
  
  const DetId& id() const { return theid; }
  const DetId& detid() const { return id(); }

  bool overflowCol() const { return overflow_(theHitColSpan); }
  
  bool overflowRow() const { return overflow_(theHitRowSpan); }
  
  bool overflow() const { return  overflowCol() || overflowRow(); }
  
  void packCol(uint16_t ymin, uint16_t yspan) {
    theMinHitCol = ymin;
    theHitColSpan = std::min(yspan, uint16_t(MAXSPAN));
  }
  void packRow(uint16_t xmin, uint16_t xspan) {
    theMinHitRow = xmin;
    theHitRowSpan = std::min(xspan, uint16_t(MAXSPAN));
  }
  
   void setClusterErrorX( float errx ) { err_x = errx; }
   void setClusterErrorY( float erry ) { err_y = erry; }
   void setClusterErrorTime( float errtime ) { err_time = errtime; }
   float getClusterErrorX() const { return err_x; }
   float getClusterErrorY() const { return err_y; }
   float getClusterErrorTime() const { return err_time; }
  
private:

   DetId theid;
   
   std::vector<uint8_t>  theHitOffset;
   std::vector<float> theHitENERGY;
   std::vector<float> theHitTIME;
   std::vector<float> theHitTIME_ERROR;
   
   uint16_t theMinHitRow=MAXPOS; // Minimum hit index in the x direction (low edge).
   uint16_t theMinHitCol=MAXPOS; // Minimum hit index in the y direction (left edge).
   uint8_t theHitRowSpan=0; // Span hit index in the x direction (low edge).
   uint8_t theHitColSpan=0; // Span hit index in the y direction (left edge).
   
   float err_x=-99999.9f;
   float err_y=-99999.9f;
   float err_time=-99999.9f;
   
   uint8_t seed_;
   
   float weighted_sum(const std::vector<float>& weights, const std::function<float (unsigned int i)>& sumFunc, const std::function<float (float,float)>& outFunc) const 
   {
     float tot=0;
     float sumW=0;
     for (unsigned int i=0; i<weights.size(); ++i)
       {
     tot += sumFunc(i);
     sumW += weights[i];
       }
     return outFunc(tot,sumW);
   }
   
   float weighted_mean(const std::vector<float>& weights, const std::function<float (unsigned int)>& value) const
   {
     auto sumFunc=[&weights,value](unsigned int i) { return weights[i]*value(i); } ;
     auto outFunc=[](float x,float y) { if (y>0) return (float)x/y; else return -999.f; };
     return weighted_sum(weights,sumFunc,outFunc);
   }
   
   float weighted_mean_error(const std::vector<float>& weights, const std::function<float (unsigned int)>& err) const
   {
     auto sumFunc=[&weights,err](unsigned int i) { return weights[i]*weights[i]*err(i)*err(i); } ;
     auto outFunc=[](float x,float y) { if (y>0) return (float)sqrt(x)/y; else return -999.f; };
     return weighted_sum(weights,sumFunc,outFunc);
   }
   
   static int overflow_(uint16_t span) { return span==uint16_t(MAXSPAN);} 
   
};


// Comparison operators  (needed by DetSetVector & SortedCollection )
inline bool operator<( const FTLCluster& one, const FTLCluster& other) {
  if(one.detid() == other.detid()) { 
    if ( one.minHitRow() < other.minHitRow() ) {
      return true;
    } else if ( one.minHitRow() > other.minHitRow() ) {
      return false;
    } else if ( one.minHitCol() < other.minHitCol() ) {
      return true;
    } else {
      return false;
    }
  }
  return one.detid() < other.detid();
}

inline bool operator<( const FTLCluster& one, const uint32_t& detid) {
  return one.detid() < detid;}

inline bool operator<( const uint32_t& detid, const FTLCluster& other) { 
  return detid < other.detid();}

#endif