From 997e9c7b6b664f52d7dade2e4c8ea46008a54f7d Mon Sep 17 00:00:00 2001 From: Fabio Cossutti Date: Tue, 8 Nov 2022 10:33:52 +0100 Subject: [PATCH 1/6] Move cluster error to weighted mean position error, based on cell size --- DataFormats/FTLRecHit/interface/FTLCluster.h | 12 ++++++++++++ .../FTLClusterizer/interface/MTDCPEBase.h | 2 ++ RecoLocalFastTime/FTLClusterizer/src/MTDCPEBase.cc | 7 ++++--- 3 files changed, 18 insertions(+), 3 deletions(-) diff --git a/DataFormats/FTLRecHit/interface/FTLCluster.h b/DataFormats/FTLRecHit/interface/FTLCluster.h index 2a45a93fee89a..51e0e26c0f9f6 100644 --- a/DataFormats/FTLRecHit/interface/FTLCluster.h +++ b/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); diff --git a/RecoLocalFastTime/FTLClusterizer/interface/MTDCPEBase.h b/RecoLocalFastTime/FTLClusterizer/interface/MTDCPEBase.h index 4330904c84224..125abb8fa5779 100644 --- a/RecoLocalFastTime/FTLClusterizer/interface/MTDCPEBase.h +++ b/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 = 1.f / std::sqrt(12.f); + protected: //--------------------------------------------------------------------------- // Data members diff --git a/RecoLocalFastTime/FTLClusterizer/src/MTDCPEBase.cc b/RecoLocalFastTime/FTLClusterizer/src/MTDCPEBase.cc index 85d7aec9a379c..8f0c6b63af812 100644 --- a/RecoLocalFastTime/FTLClusterizer/src/MTDCPEBase.cc +++ b/RecoLocalFastTime/FTLClusterizer/src/MTDCPEBase.cc @@ -66,10 +66,11 @@ LocalPoint MTDCPEBase::localPosition(DetParam const& dp, ClusterParam& cp) const } 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); + return dp.theTopol->localError(pos, posErr); } MTDCPEBase::TimeValue MTDCPEBase::clusterTime(DetParam const& dp, ClusterParam& cp) const { From 049e697688cef832ce5769844f607aaf1b83d6d0 Mon Sep 17 00:00:00 2001 From: Fabio Cossutti Date: Tue, 8 Nov 2022 15:55:05 +0100 Subject: [PATCH 2/6] Fix name for left/right sample according to DIGI convention --- .../plugins/BTLUncalibRecHitAlgo.cc | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/RecoLocalFastTime/FTLCommonAlgos/plugins/BTLUncalibRecHitAlgo.cc b/RecoLocalFastTime/FTLCommonAlgos/plugins/BTLUncalibRecHitAlgo.cc index dbf10474d9e4d..7ef8da92c0686 100644 --- a/RecoLocalFastTime/FTLCommonAlgos/plugins/BTLUncalibRecHitAlgo.cc +++ b/RecoLocalFastTime/FTLCommonAlgos/plugins/BTLUncalibRecHitAlgo.cc @@ -53,14 +53,14 @@ FTLUncalibratedRecHit BTLUncalibRecHitAlgo::makeRecHit(const BTLDataFrame& dataF unsigned char flag = 0; - const auto& sampleLeft = dataFrame.sample(0); - const auto& sampleRight = dataFrame.sample(1); + const auto& sampleRight = dataFrame.sample(0); + const auto& sampleLeft = dataFrame.sample(1); double nHits = 0.; - if (sampleLeft.data() > 0) { - amplitude.first = float(sampleLeft.data()) * adcLSB_; - time.first = float(sampleLeft.toa()) * toaLSBToNS_; + if (sampleRight.data() > 0) { + amplitude.first = float(sampleRight.data()) * adcLSB_; + time.first = float(sampleRight.toa()) * toaLSBToNS_; nHits += 1.; @@ -70,9 +70,9 @@ FTLUncalibratedRecHit BTLUncalibRecHitAlgo::makeRecHit(const BTLDataFrame& dataF } // --- If available, reconstruct the amplitude and time of the second SiPM - if (sampleRight.data() > 0) { - amplitude.second = sampleRight.data() * adcLSB_; - time.second = sampleRight.toa() * toaLSBToNS_; + if (sampleLeft.data() > 0) { + amplitude.second = float(sampleLeft.data()) * adcLSB_; + time.second = float(sampleLeft.toa()) * toaLSBToNS_; nHits += 1.; @@ -94,10 +94,10 @@ FTLUncalibratedRecHit BTLUncalibRecHitAlgo::makeRecHit(const BTLDataFrame& dataF float positionError = BTLRecHitsErrorEstimatorIM::positionError(); LogDebug("BTLUncalibRecHit") << "ADC+: set the charge to: (" << amplitude.first << ", " << amplitude.second << ") (" - << sampleLeft.data() << ", " << sampleRight.data() << " " << adcLSB_ << ' ' + << sampleRight.data() << ", " << sampleLeft.data() << " " << adcLSB_ << ' ' << std::endl; LogDebug("BTLUncalibRecHit") << "TDC+: set the time to: (" << time.first << ", " << time.second << ") (" - << sampleLeft.toa() << ", " << sampleRight.toa() << " " << toaLSBToNS_ << ' ' + << sampleRight.toa() << ", " << sampleLeft.toa() << " " << toaLSBToNS_ << ' ' << std::endl; return FTLUncalibratedRecHit( From 381ae972e203c9a3cbe2a69471417d9fea75f1b6 Mon Sep 17 00:00:00 2001 From: Fabio Cossutti Date: Thu, 10 Nov 2022 14:49:56 +0100 Subject: [PATCH 3/6] Propagate detailed position in BTL crystal to clusters --- DataFormats/FTLRecHit/interface/FTLCluster.h | 8 +-- DataFormats/FTLRecHit/src/classes_def.xml | 3 +- .../FTLClusterizer/interface/MTDArrayBuffer.h | 32 ++++++++--- .../FTLClusterizer/src/MTDCPEBase.cc | 20 ++++++- .../src/MTDThresholdClusterizer.cc | 54 ++++++++++++++----- 5 files changed, 92 insertions(+), 25 deletions(-) diff --git a/DataFormats/FTLRecHit/interface/FTLCluster.h b/DataFormats/FTLRecHit/interface/FTLCluster.h index 51e0e26c0f9f6..473a856dd1498 100644 --- a/DataFormats/FTLRecHit/interface/FTLCluster.h +++ b/DataFormats/FTLRecHit/interface/FTLCluster.h @@ -205,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: @@ -225,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_; diff --git a/DataFormats/FTLRecHit/src/classes_def.xml b/DataFormats/FTLRecHit/src/classes_def.xml index 8288abbe58e6f..3e786d54f51f4 100644 --- a/DataFormats/FTLRecHit/src/classes_def.xml +++ b/DataFormats/FTLRecHit/src/classes_def.xml @@ -23,7 +23,8 @@ - + + diff --git a/RecoLocalFastTime/FTLClusterizer/interface/MTDArrayBuffer.h b/RecoLocalFastTime/FTLClusterizer/interface/MTDArrayBuffer.h index 331dab2451ea1..5fc57a0fadc3d 100644 --- a/RecoLocalFastTime/FTLClusterizer/interface/MTDArrayBuffer.h +++ b/RecoLocalFastTime/FTLClusterizer/interface/MTDArrayBuffer.h @@ -41,6 +41,9 @@ class MTDArrayBuffer { inline GlobalPoint global_point(uint row, uint col) const; inline GlobalPoint global_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; } @@ -55,6 +58,7 @@ class MTDArrayBuffer { set_time_error(row, col, 0.); set_local_error(row, col, le_n); set_global_point(row, col, gp_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 GlobalPoint& global_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 GlobalPoint& global_point, + float xpos); inline void set_subDet(uint row, uint col, GeomDetEnumerators::Location subDet); inline void set_subDet(const FTLCluster::FTLHitPos&, GeomDetEnumerators::Location subDet); @@ -93,6 +99,9 @@ class MTDArrayBuffer { 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. @@ -106,6 +115,7 @@ class MTDArrayBuffer { std::vector hitTimeError_vec; std::vector hitGP_vec; std::vector hitLE_vec; + std::vector 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; } @@ -153,6 +164,9 @@ LocalError MTDArrayBuffer::local_error(const FTLCluster::FTLHitPos& pix) const { 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)]; } +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, GeomDetEnumerators::Location subDet, @@ -160,13 +174,15 @@ void MTDArrayBuffer::set(uint row, float time, float time_error, const LocalError& local_error, - const GlobalPoint& global_point) { + const GlobalPoint& global_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; hitLE_vec[index(row, col)] = local_error; + xpos_vec[index(row, col)] = xpos; } void MTDArrayBuffer::set(const FTLCluster::FTLHitPos& pix, GeomDetEnumerators::Location subDet, @@ -174,8 +190,9 @@ void MTDArrayBuffer::set(const FTLCluster::FTLHitPos& pix, 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 GlobalPoint& global_point, + float xpos) { + set(pix.row(), pix.col(), subDet, energy, time, time_error, local_error, global_point, xpos); } void MTDArrayBuffer::set_subDet(uint row, uint col, GeomDetEnumerators::Location subDet) { @@ -209,4 +226,7 @@ void MTDArrayBuffer::set_local_error(const FTLCluster::FTLHitPos& pix, const Loc 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 diff --git a/RecoLocalFastTime/FTLClusterizer/src/MTDCPEBase.cc b/RecoLocalFastTime/FTLClusterizer/src/MTDCPEBase.cc index 8f0c6b63af812..4f44f56b18297 100644 --- a/RecoLocalFastTime/FTLClusterizer/src/MTDCPEBase.cc +++ b/RecoLocalFastTime/FTLClusterizer/src/MTDCPEBase.cc @@ -62,7 +62,13 @@ 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 { @@ -70,7 +76,17 @@ LocalError MTDCPEBase::localError(DetParam const& dp, ClusterParam& cp) const { float sigma2 = cp.theCluster->positionError(sigma_flat); sigma2 *= sigma2; MeasurementError posErr(sigma2, 0, sigma2); - return dp.theTopol->localError(pos, posErr); + 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 { diff --git a/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc b/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc index f27ea64d6ccd7..dc46838077b1e 100644 --- a/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc +++ b/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc @@ -198,6 +198,7 @@ void MTDThresholdClusterizer::copy_to_buffer(RecHitIterator itr, const MTDGeomet float time = itr->time(); float timeError = itr->timeError(); float position = itr->position(); + float xpos = 0.f; // position is the longitudinal offset that should be added into local x for bars in phi geometry LocalError local_error(0, 0, 0); GlobalPoint global_point(0, 0, 0); @@ -208,14 +209,13 @@ void MTDThresholdClusterizer::copy_to_buffer(RecHitIterator itr, const MTDGeomet const auto& det = geom->idToDet(geoId); const ProxyMTDTopology& topoproxy = static_cast(det->topology()); const RectangularMTDTopology& topol = static_cast(topoproxy.specificTopology()); - MeasurementPoint mp(row, col); - LocalPoint lp_ctr = topol.localPosition(mp); - LocalPoint lp(lp_ctr.x() + position + topol.pitch().first * 0.5f, lp_ctr.y(), lp_ctr.z()); - // local coordinates of BTL module locates RecHits on the left edge of the bar (-9.2, -3.067, 3.067) - // (position + topol.pitch().first/2.0) is the distance from the left edge to the Hit point - global_point = det->toGlobal(lp); - BTLRecHitsErrorEstimatorIM btlError(det, lp); + + LocalPoint lp_pixel(position, 0, 0); + LocalPoint lp_module = topol.pixelToModuleLocalPoint(lp_pixel, row, col); + global_point = det->toGlobal(lp_module); + BTLRecHitsErrorEstimatorIM btlError(det, lp_module); local_error = btlError.localError(); + xpos = lp_module.x(); } else if (mtdId.mtdSubDetector() == MTDDetId::ETL) { subDet = GeomDetEnumerators::endcap; ETLDetId id = itr->id(); @@ -224,14 +224,18 @@ void MTDThresholdClusterizer::copy_to_buffer(RecHitIterator itr, const MTDGeomet const ProxyMTDTopology& topoproxy = static_cast(det->topology()); const RectangularMTDTopology& topol = static_cast(topoproxy.specificTopology()); - MeasurementPoint mp(row, col); - LocalPoint lp = topol.localPosition(mp); - global_point = det->toGlobal(lp); + LocalPoint lp_pixel(0, 0, 0); + LocalPoint lp_module = topol.pixelToModuleLocalPoint(lp_pixel, row, col); + global_point = det->toGlobal(lp_module); } - LogDebug("MTDThresholdClusterizer") << "ROW " << row << " COL " << col << " ENERGY " << energy << " TIME " << time; + LogDebug("MTDThresholdClusterizer") << "DetId " << mtdId.rawId() << " subd " << mtdId.mtdSubDetector() << " row/col " + << row << " / " << col << " energy " << energy << " time " << time + << " time error " << timeError << " global_point " << global_point.x() << " " + << global_point.y() << " " << global_point.z() << " local error " + << local_error.xx() << " " << local_error.yy() << " xpos " << xpos; if (energy > theHitThreshold) { - theBuffer.set(row, col, subDet, energy, time, timeError, local_error, global_point); + theBuffer.set(row, col, subDet, energy, time, timeError, local_error, global_point, xpos); if (energy > theSeedThreshold) theSeeds.push_back(FTLCluster::FTLHitPos(row, col)); //sort seeds? @@ -256,6 +260,14 @@ FTLCluster MTDThresholdClusterizer::make_cluster(const FTLCluster::FTLHitPos& hi AccretionCluster acluster; acluster.add(hit, seed_energy, seed_time, seed_time_error); + // for BTL position along crystals add auxiliary vectors + std::array pixel_x{{-99999.}}; + std::array pixel_errx2{{-99999.}}; + if (seed_subdet == GeomDetEnumerators::barrel) { + pixel_x[acluster.top()] = theBuffer.xpos(hit.row(), hit.col()); + pixel_errx2[acluster.top()] = seed_error_xx; + } + bool stopClus = false; //Here we search all hits adjacent to all hits in the cluster. while (!acluster.empty() && !stopClus) { @@ -288,6 +300,10 @@ FTLCluster MTDThresholdClusterizer::make_cluster(const FTLCluster::FTLHitPos& hi stopClus = true; break; } + if (theBuffer.subDet(r, c) == GeomDetEnumerators::barrel) { + pixel_x[curInd] = theBuffer.xpos(r, c); + pixel_errx2[curInd] = theBuffer.local_error(r, c).xx(); + } theBuffer.clear(newhit); } } @@ -303,5 +319,19 @@ FTLCluster MTDThresholdClusterizer::make_cluster(const FTLCluster::FTLHitPos& hi acluster.y.data(), acluster.xmin, acluster.ymin); + + // For BTL compute the optimal position along crystal and uncertainty on it in absolute length units + + if (seed_subdet == GeomDetEnumerators::barrel) { + float sumW(0.f), sumXW(0.f), sumXW2(0.f); + for (unsigned int index = 0; index < acluster.top(); index++) { + sumW += acluster.energy[index]; + sumXW += acluster.energy[index] * pixel_x[index]; + sumXW2 += acluster.energy[index] * acluster.energy[index] * pixel_errx2[index]; + } + cluster.setClusterPosX(sumXW / sumW); + cluster.setClusterErrorX(std::sqrt(sumXW2 / sumW / sumW)); + } + return cluster; } From eb2e7ee8bb2d12ed911299153447d1843c18fc0b Mon Sep 17 00:00:00 2001 From: Fabio Cossutti Date: Fri, 18 Nov 2022 16:49:11 +0100 Subject: [PATCH 4/6] Fix x position clustering, move GlobaPoint to LocalPoint, update debugging statements --- .../FTLClusterizer/interface/MTDArrayBuffer.h | 38 +++++++++---------- .../src/MTDThresholdClusterizer.cc | 36 +++++++++--------- 2 files changed, 38 insertions(+), 36 deletions(-) diff --git a/RecoLocalFastTime/FTLClusterizer/interface/MTDArrayBuffer.h b/RecoLocalFastTime/FTLClusterizer/interface/MTDArrayBuffer.h index 5fc57a0fadc3d..0a67b7cdffcf7 100644 --- a/RecoLocalFastTime/FTLClusterizer/interface/MTDArrayBuffer.h +++ b/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 @@ -38,8 +38,8 @@ 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; @@ -51,13 +51,13 @@ 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()); } @@ -69,7 +69,7 @@ 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, @@ -77,7 +77,7 @@ class MTDArrayBuffer { 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); @@ -93,8 +93,8 @@ 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); @@ -113,7 +113,7 @@ class MTDArrayBuffer { std::vector hitEnergy_vec; std::vector hitTime_vec; std::vector hitTimeError_vec; - std::vector hitGP_vec; + std::vector hitGP_vec; std::vector hitLE_vec; std::vector xpos_vec; uint nrows; @@ -161,8 +161,8 @@ 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)]; } @@ -174,13 +174,13 @@ 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; } @@ -190,9 +190,9 @@ void MTDArrayBuffer::set(const FTLCluster::FTLHitPos& pix, float time, float time_error, const LocalError& local_error, - const GlobalPoint& global_point, + const LocalPoint& local_point, float xpos) { - set(pix.row(), pix.col(), subDet, energy, time, time_error, local_error, global_point, 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) { @@ -216,9 +216,9 @@ 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; } diff --git a/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc b/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc index dc46838077b1e..ed621c8592d71 100644 --- a/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc +++ b/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc @@ -161,7 +161,9 @@ void MTDThresholdClusterizer::clusterize(const FTLRecHitCollection& input, if (cluster.energy() > theClusterThreshold) { LogDebug("MTDThresholdClusterizer") << "putting in this cluster " << i << " #hits:" << cluster.size() << " E:" << cluster.energy() - << " T:" << cluster.time() << " X:" << cluster.x() << " Y:" << cluster.y(); + << " T +/- DT:" << cluster.time() << " +/- " << cluster.timeError() << " X:" << cluster.x() + << " Y:" << cluster.y() << " xpos +/- err " << cluster.getClusterPosX() << " +/- " + << cluster.getClusterErrorX(); clustersOnDet.push_back(cluster); } } @@ -200,8 +202,8 @@ void MTDThresholdClusterizer::copy_to_buffer(RecHitIterator itr, const MTDGeomet float position = itr->position(); float xpos = 0.f; // position is the longitudinal offset that should be added into local x for bars in phi geometry + LocalPoint local_point(0, 0, 0); LocalError local_error(0, 0, 0); - GlobalPoint global_point(0, 0, 0); if (mtdId.mtdSubDetector() == MTDDetId::BTL) { subDet = GeomDetEnumerators::barrel; BTLDetId id = itr->id(); @@ -211,11 +213,10 @@ void MTDThresholdClusterizer::copy_to_buffer(RecHitIterator itr, const MTDGeomet const RectangularMTDTopology& topol = static_cast(topoproxy.specificTopology()); LocalPoint lp_pixel(position, 0, 0); - LocalPoint lp_module = topol.pixelToModuleLocalPoint(lp_pixel, row, col); - global_point = det->toGlobal(lp_module); - BTLRecHitsErrorEstimatorIM btlError(det, lp_module); + local_point = topol.pixelToModuleLocalPoint(lp_pixel, row, col); + BTLRecHitsErrorEstimatorIM btlError(det, local_point); local_error = btlError.localError(); - xpos = lp_module.x(); + xpos = local_point.x(); } else if (mtdId.mtdSubDetector() == MTDDetId::ETL) { subDet = GeomDetEnumerators::endcap; ETLDetId id = itr->id(); @@ -225,17 +226,17 @@ void MTDThresholdClusterizer::copy_to_buffer(RecHitIterator itr, const MTDGeomet const RectangularMTDTopology& topol = static_cast(topoproxy.specificTopology()); LocalPoint lp_pixel(0, 0, 0); - LocalPoint lp_module = topol.pixelToModuleLocalPoint(lp_pixel, row, col); - global_point = det->toGlobal(lp_module); + local_point = topol.pixelToModuleLocalPoint(lp_pixel, row, col); } LogDebug("MTDThresholdClusterizer") << "DetId " << mtdId.rawId() << " subd " << mtdId.mtdSubDetector() << " row/col " << row << " / " << col << " energy " << energy << " time " << time - << " time error " << timeError << " global_point " << global_point.x() << " " - << global_point.y() << " " << global_point.z() << " local error " - << local_error.xx() << " " << local_error.yy() << " xpos " << xpos; + << " time error " << timeError << " local_point " << local_point.x() << " " + << local_point.y() << " " << local_point.z() << " local error " + << std::sqrt(local_error.xx()) << " " << std::sqrt(local_error.yy()) << " xpos " + << xpos; if (energy > theHitThreshold) { - theBuffer.set(row, col, subDet, energy, time, timeError, local_error, global_point, xpos); + theBuffer.set(row, col, subDet, energy, time, timeError, local_error, local_point, xpos); if (energy > theSeedThreshold) theSeeds.push_back(FTLCluster::FTLHitPos(row, col)); //sort seeds? @@ -252,9 +253,10 @@ FTLCluster MTDThresholdClusterizer::make_cluster(const FTLCluster::FTLHitPos& hi float seed_energy = theBuffer.energy(hit.row(), hit.col()); float seed_time = theBuffer.time(hit.row(), hit.col()); float seed_time_error = theBuffer.time_error(hit.row(), hit.col()); - auto const seedPoint = theBuffer.global_point(hit.row(), hit.col()); + auto const seedPoint = theBuffer.local_point(hit.row(), hit.col()); double seed_error_xx = theBuffer.local_error(hit.row(), hit.col()).xx(); double seed_error_yy = theBuffer.local_error(hit.row(), hit.col()).yy(); + double seed_xpos = theBuffer.xpos(hit.row(), hit.col()); theBuffer.clear(hit); AccretionCluster acluster; @@ -264,7 +266,7 @@ FTLCluster MTDThresholdClusterizer::make_cluster(const FTLCluster::FTLHitPos& hi std::array pixel_x{{-99999.}}; std::array pixel_errx2{{-99999.}}; if (seed_subdet == GeomDetEnumerators::barrel) { - pixel_x[acluster.top()] = theBuffer.xpos(hit.row(), hit.col()); + pixel_x[acluster.top()] = seed_xpos; pixel_errx2[acluster.top()] = seed_error_xx; } @@ -289,7 +291,7 @@ FTLCluster MTDThresholdClusterizer::make_cluster(const FTLCluster::FTLHitPos& hi double hit_error_xx = theBuffer.local_error(r, c).xx(); double hit_error_yy = theBuffer.local_error(r, c).yy(); if (thePositionThreshold > 0) { - if (((theBuffer.global_point(r, c) - seedPoint).mag2()) > + if (((theBuffer.local_point(r, c) - seedPoint).mag2()) > thePositionThreshold * thePositionThreshold * (hit_error_xx + seed_error_xx + hit_error_yy + seed_error_yy)) continue; @@ -301,8 +303,8 @@ FTLCluster MTDThresholdClusterizer::make_cluster(const FTLCluster::FTLHitPos& hi break; } if (theBuffer.subDet(r, c) == GeomDetEnumerators::barrel) { - pixel_x[curInd] = theBuffer.xpos(r, c); - pixel_errx2[curInd] = theBuffer.local_error(r, c).xx(); + pixel_x[acluster.top()] = theBuffer.xpos(r, c); + pixel_errx2[acluster.top()] = theBuffer.local_error(r, c).xx(); } theBuffer.clear(newhit); } From 812db3e924fc92865f7c39ee4551114c6d8d42f2 Mon Sep 17 00:00:00 2001 From: Fabio Cossutti Date: Tue, 22 Nov 2022 10:01:32 +0100 Subject: [PATCH 5/6] Minor technical adjustments --- .../FTLClusterizer/interface/MTDCPEBase.h | 2 +- .../src/MTDThresholdClusterizer.cc | 19 +++++++++++-------- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/RecoLocalFastTime/FTLClusterizer/interface/MTDCPEBase.h b/RecoLocalFastTime/FTLClusterizer/interface/MTDCPEBase.h index 125abb8fa5779..44f38ecb84ace 100644 --- a/RecoLocalFastTime/FTLClusterizer/interface/MTDCPEBase.h +++ b/RecoLocalFastTime/FTLClusterizer/interface/MTDCPEBase.h @@ -80,7 +80,7 @@ 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 = 1.f / std::sqrt(12.f); + static constexpr float sigma_flat = 0.2886751f; // 1.f / std::sqrt(12.f); protected: //--------------------------------------------------------------------------- diff --git a/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc b/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc index ed621c8592d71..a120e01e30a14 100644 --- a/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc +++ b/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc @@ -74,7 +74,7 @@ bool MTDThresholdClusterizer::setup(const MTDGeometry* geom, const MTDTopology* if (nrows > theBuffer.rows() || ncols > theBuffer.columns()) { // change only when a larger is needed // Resize the buffer - theBuffer.setSize(nrows + 1, ncols + 1); // +1 needed for MTD + theBuffer.setSize(nrows, ncols); bufferAlreadySet = true; } @@ -282,20 +282,23 @@ FTLCluster MTDThresholdClusterizer::make_cluster(const FTLCluster::FTLHitPos& hi for (auto r = std::max(0, int(acluster.x[curInd] - 1)); r < std::min(int(acluster.x[curInd] + 2), int(theBuffer.rows())) && !stopClus; ++r) { + LogDebug("MTDThresholdClusterizer") + << "Clustering subdet " << seed_subdet << " around " << curInd << " X,Y = " << acluster.x[curInd] << "," + << acluster.y[curInd] << " r,c = " << r << "," << c << " energy,time = " << theBuffer.energy(r, c) << " " + << theBuffer.time(r, c); if (theBuffer.energy(r, c) > theHitThreshold) { if (std::abs(theBuffer.time(r, c) - seed_time) > theTimeThreshold * sqrt(theBuffer.time_error(r, c) * theBuffer.time_error(r, c) + seed_time_error * seed_time_error)) continue; - if ((seed_subdet == GeomDetEnumerators::barrel) && (theBuffer.subDet(r, c) == GeomDetEnumerators::barrel)) { + if ((seed_subdet == GeomDetEnumerators::barrel) && (theBuffer.subDet(r, c) == GeomDetEnumerators::barrel) && + (thePositionThreshold > 0)) { double hit_error_xx = theBuffer.local_error(r, c).xx(); double hit_error_yy = theBuffer.local_error(r, c).yy(); - if (thePositionThreshold > 0) { - if (((theBuffer.local_point(r, c) - seedPoint).mag2()) > - thePositionThreshold * thePositionThreshold * - (hit_error_xx + seed_error_xx + hit_error_yy + seed_error_yy)) - continue; - } + if (((theBuffer.local_point(r, c) - seedPoint).mag2()) > + thePositionThreshold * thePositionThreshold * + (hit_error_xx + seed_error_xx + hit_error_yy + seed_error_yy)) + continue; } FTLCluster::FTLHitPos newhit(r, c); if (!acluster.add(newhit, theBuffer.energy(r, c), theBuffer.time(r, c), theBuffer.time_error(r, c))) { From b314558a28e1516f95af6420e0628cc8cd9f59cc Mon Sep 17 00:00:00 2001 From: Fabio Cossutti Date: Sun, 27 Nov 2022 18:04:51 +0100 Subject: [PATCH 6/6] Comment by M. Nguyen --- RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc b/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc index a120e01e30a14..9e20fe8c0782d 100644 --- a/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc +++ b/RecoLocalFastTime/FTLClusterizer/src/MTDThresholdClusterizer.cc @@ -335,7 +335,7 @@ FTLCluster MTDThresholdClusterizer::make_cluster(const FTLCluster::FTLHitPos& hi sumXW2 += acluster.energy[index] * acluster.energy[index] * pixel_errx2[index]; } cluster.setClusterPosX(sumXW / sumW); - cluster.setClusterErrorX(std::sqrt(sumXW2 / sumW / sumW)); + cluster.setClusterErrorX(std::sqrt(sumXW2) / sumW); } return cluster;