MTD reconstruction: fix local cluster position and uncertainty determination

DataFormats/FTLRecHit/interface/FTLCluster.h

@@ -127,6 +127,18 @@ class FTLCluster {
     return weighted_mean(this->theHitENERGY, y_pos);
   }
 
+  inline float positionError(const float sigmaPos) const {
+    float sumW2(0.f), sumW(0.f);
+    for (const auto& hitW : theHitENERGY) {
+      sumW2 += hitW * hitW;
+      sumW += hitW;
+    }
+    if (sumW > 0)
+      return sigmaPos * std::sqrt(sumW2) / sumW;
+    else
+      return -999.f;
+  }
+
   inline float time() const {
     auto t = [this](unsigned int i) { return this->theHitTIME[i]; };
     return weighted_mean(this->theHitENERGY, t);
@@ -193,11 +205,11 @@ class FTLCluster {
     theHitRowSpan = std::min(xspan, uint16_t(MAXSPAN));
   }
 
+  void setClusterPosX(float posx) { pos_x = posx; }
   void setClusterErrorX(float errx) { err_x = errx; }
-  void setClusterErrorY(float erry) { err_y = erry; }
   void setClusterErrorTime(float errtime) { err_time = errtime; }
+  float getClusterPosX() const { return pos_x; }
   float getClusterErrorX() const { return err_x; }
-  float getClusterErrorY() const { return err_y; }
   float getClusterErrorTime() const { return err_time; }
 
 private:
@@ -213,8 +225,8 @@ class FTLCluster {
   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 pos_x = -99999.9f;  // For pixels with internal position information in one coordinate (i.e. BTL crystals)
+  float err_x = -99999.9f;  // For pixels with internal position information in one coordinate (i.e. BTL crystals)
   float err_time = -99999.9f;
 
   uint8_t seed_;

DataFormats/FTLRecHit/src/classes_def.xml

@@ -23,7 +23,8 @@
   <class name="edm::RefVector<FTLRecHitCollection>"/>
   <class name="edm::RefProd<FTLRecHitCollection>"/>
 
-  <class name="FTLCluster" ClassVersion="3">
+  <class name="FTLCluster" ClassVersion="4">
+   <version ClassVersion="4" checksum="291411646"/>
     <version ClassVersion="3" checksum="2080981492"/>
   </class>
   <class name="std::vector<FTLCluster>"/>

RecoLocalFastTime/FTLClusterizer/interface/MTDArrayBuffer.h

@@ -11,7 +11,7 @@
 #include "DataFormats/FTLRecHit/interface/FTLCluster.h"
 
 #include "DataFormats/GeometrySurface/interface/LocalError.h"
-#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
+#include "DataFormats/GeometryVector/interface/LocalPoint.h"
 #include "Geometry/CommonDetUnit/interface/GeomDetEnumerators.h"
 
 #include <vector>
@@ -38,8 +38,11 @@ class MTDArrayBuffer {
 
   inline LocalError local_error(uint row, uint col) const;
   inline LocalError local_error(const FTLCluster::FTLHitPos&) const;
-  inline GlobalPoint global_point(uint row, uint col) const;
-  inline GlobalPoint global_point(const FTLCluster::FTLHitPos&) const;
+  inline LocalPoint local_point(uint row, uint col) const;
+  inline LocalPoint local_point(const FTLCluster::FTLHitPos&) const;
+
+  inline float xpos(uint row, uint col) const;
+  inline float xpos(const FTLCluster::FTLHitPos&) const;
 
   inline uint rows() const { return nrows; }
   inline uint columns() const { return ncols; }
@@ -48,13 +51,14 @@ class MTDArrayBuffer {
 
   inline void clear(uint row, uint col) {
     LocalError le_n(0, 0, 0);
-    GlobalPoint gp_n(0, 0, 0);
+    LocalPoint lp_n(0, 0, 0);
     set_subDet(row, col, GeomDetEnumerators::invalidLoc);
     set_energy(row, col, 0.);
     set_time(row, col, 0.);
     set_time_error(row, col, 0.);
     set_local_error(row, col, le_n);
-    set_global_point(row, col, gp_n);
+    set_local_point(row, col, lp_n);
+    set_xpos(row, col, 0.);
   }
   inline void clear(const FTLCluster::FTLHitPos& pos) { clear(pos.row(), pos.col()); }
 
@@ -65,14 +69,16 @@ class MTDArrayBuffer {
                   float time,
                   float time_error,
                   const LocalError& local_error,
-                  const GlobalPoint& global_point);
+                  const LocalPoint& local_point,
+                  float xpos);
   inline void set(const FTLCluster::FTLHitPos&,
                   GeomDetEnumerators::Location subDet,
                   float energy,
                   float time,
                   float time_error,
                   const LocalError& local_error,
-                  const GlobalPoint& global_point);
+                  const LocalPoint& local_point,
+                  float xpos);
 
   inline void set_subDet(uint row, uint col, GeomDetEnumerators::Location subDet);
   inline void set_subDet(const FTLCluster::FTLHitPos&, GeomDetEnumerators::Location subDet);
@@ -87,12 +93,15 @@ class MTDArrayBuffer {
   inline void set_time_error(uint row, uint col, float time_error);
   inline void set_time_error(const FTLCluster::FTLHitPos&, float time_error);
 
-  inline void set_global_point(uint row, uint col, const GlobalPoint& gp);
-  inline void set_global_point(const FTLCluster::FTLHitPos&, const GlobalPoint& gp);
+  inline void set_local_point(uint row, uint col, const LocalPoint& lp);
+  inline void set_local_point(const FTLCluster::FTLHitPos&, const LocalPoint& lp);
 
   inline void set_local_error(uint row, uint col, const LocalError& le);
   inline void set_local_error(const FTLCluster::FTLHitPos&, const LocalError& le);
 
+  inline void set_xpos(uint row, uint col, float xpos);
+  inline void set_xpos(const FTLCluster::FTLHitPos&, float xpos);
+
   uint size() const { return hitEnergy_vec.size(); }
 
   /// Definition of indexing within the buffer.
@@ -104,8 +113,9 @@ class MTDArrayBuffer {
   std::vector<float> hitEnergy_vec;
   std::vector<float> hitTime_vec;
   std::vector<float> hitTimeError_vec;
-  std::vector<GlobalPoint> hitGP_vec;
+  std::vector<LocalPoint> hitGP_vec;
   std::vector<LocalError> hitLE_vec;
+  std::vector<float> xpos_vec;
   uint nrows;
   uint ncols;
 };
@@ -117,6 +127,7 @@ MTDArrayBuffer::MTDArrayBuffer(uint rows, uint cols)
       hitTimeError_vec(rows * cols, 0),
       hitGP_vec(rows * cols),
       hitLE_vec(rows * cols),
+      xpos_vec(rows * cols),
       nrows(rows),
       ncols(cols) {}
 
@@ -127,7 +138,7 @@ void MTDArrayBuffer::setSize(uint rows, uint cols) {
   hitTimeError_vec.resize(rows * cols, 0);
   hitGP_vec.resize(rows * cols);
   hitLE_vec.resize(rows * cols);
-  nrows = rows;
+  xpos_vec.resize(rows * cols), nrows = rows;
   ncols = cols;
 }
 
@@ -150,8 +161,11 @@ float MTDArrayBuffer::time_error(const FTLCluster::FTLHitPos& pix) const { retur
 LocalError MTDArrayBuffer::local_error(uint row, uint col) const { return hitLE_vec[index(row, col)]; }
 LocalError MTDArrayBuffer::local_error(const FTLCluster::FTLHitPos& pix) const { return hitLE_vec[index(pix)]; }
 
-GlobalPoint MTDArrayBuffer::global_point(uint row, uint col) const { return hitGP_vec[index(row, col)]; }
-GlobalPoint MTDArrayBuffer::global_point(const FTLCluster::FTLHitPos& pix) const { return hitGP_vec[index(pix)]; }
+LocalPoint MTDArrayBuffer::local_point(uint row, uint col) const { return hitGP_vec[index(row, col)]; }
+LocalPoint MTDArrayBuffer::local_point(const FTLCluster::FTLHitPos& pix) const { return hitGP_vec[index(pix)]; }
+
+float MTDArrayBuffer::xpos(uint row, uint col) const { return xpos_vec[index(row, col)]; }
+float MTDArrayBuffer::xpos(const FTLCluster::FTLHitPos& pix) const { return xpos_vec[index(pix)]; }
 
 void MTDArrayBuffer::set(uint row,
                          uint col,
@@ -160,22 +174,25 @@ void MTDArrayBuffer::set(uint row,
                          float time,
                          float time_error,
                          const LocalError& local_error,
-                         const GlobalPoint& global_point) {
+                         const LocalPoint& local_point,
+                         float xpos) {
   hitSubDet_vec[index(row, col)] = subDet;
   hitEnergy_vec[index(row, col)] = energy;
   hitTime_vec[index(row, col)] = time;
   hitTimeError_vec[index(row, col)] = time_error;
-  hitGP_vec[index(row, col)] = global_point;
+  hitGP_vec[index(row, col)] = local_point;
   hitLE_vec[index(row, col)] = local_error;
+  xpos_vec[index(row, col)] = xpos;
 }
 void MTDArrayBuffer::set(const FTLCluster::FTLHitPos& pix,
                          GeomDetEnumerators::Location subDet,
                          float energy,
                          float time,
                          float time_error,
                          const LocalError& local_error,
-                         const GlobalPoint& global_point) {
-  set(pix.row(), pix.col(), subDet, energy, time, time_error, local_error, global_point);
+                         const LocalPoint& local_point,
+                         float xpos) {
+  set(pix.row(), pix.col(), subDet, energy, time, time_error, local_error, local_point, xpos);
 }
 
 void MTDArrayBuffer::set_subDet(uint row, uint col, GeomDetEnumerators::Location subDet) {
@@ -199,14 +216,17 @@ void MTDArrayBuffer::set_time_error(const FTLCluster::FTLHitPos& pix, float time
   hitTimeError_vec[index(pix)] = time_error;
 }
 
-void MTDArrayBuffer::set_global_point(uint row, uint col, const GlobalPoint& gp) { hitGP_vec[index(row, col)] = gp; }
-void MTDArrayBuffer::set_global_point(const FTLCluster::FTLHitPos& pix, const GlobalPoint& gp) {
-  hitGP_vec[index(pix)] = gp;
+void MTDArrayBuffer::set_local_point(uint row, uint col, const LocalPoint& lp) { hitGP_vec[index(row, col)] = lp; }
+void MTDArrayBuffer::set_local_point(const FTLCluster::FTLHitPos& pix, const LocalPoint& lp) {
+  hitGP_vec[index(pix)] = lp;
 }
 
 void MTDArrayBuffer::set_local_error(uint row, uint col, const LocalError& le) { hitLE_vec[index(row, col)] = le; }
 void MTDArrayBuffer::set_local_error(const FTLCluster::FTLHitPos& pix, const LocalError& le) {
   hitLE_vec[index(pix)] = le;
 }
 
+void MTDArrayBuffer::set_xpos(uint row, uint col, float xpos) { xpos_vec[index(row, col)] = xpos; }
+void MTDArrayBuffer::set_xpos(const FTLCluster::FTLHitPos& pix, float xpos) { xpos_vec[index(pix)] = xpos; }
+
 #endif

RecoLocalFastTime/FTLClusterizer/interface/MTDCPEBase.h

@@ -80,6 +80,8 @@ class MTDCPEBase : public MTDClusterParameterEstimator {
   virtual TimeValue clusterTime(DetParam const& dp, ClusterParam& cp) const;
   virtual TimeValueError clusterTimeError(DetParam const& dp, ClusterParam& cp) const;
 
+  static constexpr float sigma_flat = 0.2886751f;  // 1.f / std::sqrt(12.f);
+
 protected:
   //---------------------------------------------------------------------------
   //  Data members

RecoLocalFastTime/FTLClusterizer/src/MTDCPEBase.cc

@@ -62,14 +62,31 @@ MTDCPEBase::DetParam const& MTDCPEBase::detParam(const GeomDetUnit& det) const {
 LocalPoint MTDCPEBase::localPosition(DetParam const& dp, ClusterParam& cp) const {
   //remember measurement point is row(col)+0.5f
   MeasurementPoint pos(cp.theCluster->x(), cp.theCluster->y());
-  return dp.theTopol->localPosition(pos);
+
+  if (cp.theCluster->getClusterErrorX() < 0.) {
+    return dp.theTopol->localPosition(pos);
+  } else {
+    LocalPoint out(cp.theCluster->getClusterPosX(), dp.theTopol->localY(pos.y()));
+    return out;
+  }
 }
 
 LocalError MTDCPEBase::localError(DetParam const& dp, ClusterParam& cp) const {
-  constexpr double one_over_twelve = 1. / 12.;
   MeasurementPoint pos(cp.theCluster->x(), cp.theCluster->y());
-  MeasurementError simpleRect(one_over_twelve, 0, one_over_twelve);
-  return dp.theTopol->localError(pos, simpleRect);
+  float sigma2 = cp.theCluster->positionError(sigma_flat);
+  sigma2 *= sigma2;
+  MeasurementError posErr(sigma2, 0, sigma2);
+  LocalError outErr = dp.theTopol->localError(pos, posErr);
+
+  if (cp.theCluster->getClusterErrorX() < 0.) {
+    return outErr;
+  } else {
+    LocalPoint outPos(cp.theCluster->getClusterPosX(), dp.theTopol->localY(pos.y()));
+    float sigmaX2 = cp.theCluster->getClusterErrorX();
+    sigmaX2 *= sigmaX2;
+    LocalError outErrDet(sigmaX2, 0, outErr.yy());
+    return outErrDet;
+  }
 }
 
 MTDCPEBase::TimeValue MTDCPEBase::clusterTime(DetParam const& dp, ClusterParam& cp) const {