File indexing completed on 2024-04-06 12:25:56
0001
0002 #include "RecoLocalFastTime/FTLClusterizer/interface/MTDArrayBuffer.h"
0003 #include "RecoLocalFastTime/FTLClusterizer/interface/MTDThresholdClusterizer.h"
0004 #include "RecoLocalFastTime/FTLClusterizer/interface/BTLRecHitsErrorEstimatorIM.h"
0005
0006 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0007 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0008
0009 #include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
0010 #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
0011 #include "Geometry/MTDCommonData/interface/MTDTopologyMode.h"
0012
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015
0016 #include <stack>
0017 #include <vector>
0018 #include <iostream>
0019 #include <atomic>
0020 #include <cmath>
0021 using namespace std;
0022
0023
0024
0025
0026 MTDThresholdClusterizer::MTDThresholdClusterizer(edm::ParameterSet const& conf)
0027 :
0028 theHitThreshold(conf.getParameter<double>("HitThreshold")),
0029 theSeedThreshold(conf.getParameter<double>("SeedThreshold")),
0030 theClusterThreshold(conf.getParameter<double>("ClusterThreshold")),
0031 theTimeThreshold(conf.getParameter<double>("TimeThreshold")),
0032 thePositionThreshold(conf.getParameter<double>("PositionThreshold")),
0033 theNumOfRows(0),
0034 theNumOfCols(0),
0035 theCurrentId(0),
0036 theBuffer(theNumOfRows, theNumOfCols),
0037 bufferAlreadySet(false) {}
0038
0039 MTDThresholdClusterizer::~MTDThresholdClusterizer() {}
0040
0041
0042 void MTDThresholdClusterizer::fillPSetDescription(edm::ParameterSetDescription& desc) {
0043 desc.add<double>("HitThreshold", 0.);
0044 desc.add<double>("SeedThreshold", 0.);
0045 desc.add<double>("ClusterThreshold", 0.);
0046 desc.add<double>("TimeThreshold", 10.);
0047 desc.add<double>("PositionThreshold", -1.0);
0048 }
0049
0050
0051
0052
0053
0054 bool MTDThresholdClusterizer::setup(const MTDGeometry* geom, const MTDTopology* topo, const DetId& id) {
0055 theCurrentId = id;
0056
0057 const auto& thedet = geom->idToDet(id);
0058 if (thedet == nullptr) {
0059 throw cms::Exception("MTDThresholdClusterizer")
0060 << "GeographicalID: " << std::hex << id.rawId() << " is invalid!" << std::dec << std::endl;
0061 return false;
0062 }
0063 const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
0064 const RectangularMTDTopology& topol = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0065
0066
0067 unsigned int nrows = topol.nrows();
0068 unsigned int ncols = topol.ncolumns();
0069
0070 theNumOfRows = nrows;
0071 theNumOfCols = ncols;
0072
0073 LogDebug("MTDThresholdClusterizer") << "Buffer size [" << theNumOfRows + 1 << "," << theNumOfCols + 1 << "]";
0074
0075 if (nrows > theBuffer.rows() || ncols > theBuffer.columns()) {
0076
0077 theBuffer.setSize(nrows, ncols);
0078 bufferAlreadySet = true;
0079 }
0080
0081 return true;
0082 }
0083
0084
0085
0086
0087
0088
0089
0090 void MTDThresholdClusterizer::clusterize(const FTLRecHitCollection& input,
0091 const MTDGeometry* geom,
0092 const MTDTopology* topo,
0093 FTLClusterCollection& output) {
0094 FTLRecHitCollection::const_iterator begin = input.begin();
0095 FTLRecHitCollection::const_iterator end = input.end();
0096
0097
0098 if (begin == end) {
0099 edm::LogInfo("MTDThresholdClusterizer") << "No hits to clusterize";
0100 return;
0101 }
0102
0103 LogDebug("MTDThresholdClusterizer") << "Input collection " << input.size();
0104 assert(output.empty());
0105
0106 std::set<unsigned> geoIds;
0107 std::multimap<unsigned, unsigned> geoIdToIdx;
0108
0109 unsigned index = 0;
0110 for (const auto& hit : input) {
0111 MTDDetId mtdId = MTDDetId(hit.detid());
0112 if (mtdId.subDetector() != MTDDetId::FastTime) {
0113 throw cms::Exception("MTDThresholdClusterizer")
0114 << "MTDDetId: " << std::hex << mtdId.rawId() << " is invalid!" << std::dec << std::endl;
0115 }
0116
0117 if (mtdId.mtdSubDetector() == MTDDetId::BTL) {
0118 BTLDetId hitId(hit.detid());
0119
0120 DetId geoId = hitId.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topo->getMTDTopologyMode()));
0121 geoIdToIdx.emplace(geoId, index);
0122 geoIds.emplace(geoId);
0123 ++index;
0124 } else if (mtdId.mtdSubDetector() == MTDDetId::ETL) {
0125 ETLDetId hitId(hit.detid());
0126 DetId geoId = hitId.geographicalId();
0127 geoIdToIdx.emplace(geoId, index);
0128 geoIds.emplace(geoId);
0129 ++index;
0130 } else {
0131 throw cms::Exception("MTDThresholdClusterizer")
0132 << "MTDDetId: " << std::hex << mtdId.rawId() << " is invalid!" << std::dec << std::endl;
0133 }
0134 }
0135
0136
0137 for (unsigned id : geoIds) {
0138
0139 if (!setup(geom, topo, DetId(id)))
0140 return;
0141
0142 auto range = geoIdToIdx.equal_range(id);
0143 LogDebug("MTDThresholdClusterizer") << "Matching Ids for " << std::hex << id << std::dec << " ["
0144 << range.first->second << "," << range.second->second << "]";
0145
0146
0147
0148 for (auto itr = range.first; itr != range.second; ++itr) {
0149 const unsigned hitidx = itr->second;
0150 copy_to_buffer(begin + hitidx, geom, topo);
0151 }
0152
0153 FTLClusterCollection::FastFiller clustersOnDet(output, id);
0154
0155 for (unsigned int i = 0; i < theSeeds.size(); i++) {
0156 if (theBuffer.energy(theSeeds[i]) > theSeedThreshold) {
0157
0158 const FTLCluster& cluster = make_cluster(theSeeds[i]);
0159
0160
0161 if (cluster.energy() > theClusterThreshold) {
0162 LogDebug("MTDThresholdClusterizer")
0163 << "putting in this cluster " << i << " #hits:" << cluster.size() << " E:" << cluster.energy()
0164 << " T +/- DT:" << cluster.time() << " +/- " << cluster.timeError() << " X:" << cluster.x()
0165 << " Y:" << cluster.y() << " xpos +/- err " << cluster.getClusterPosX() << " +/- "
0166 << cluster.getClusterErrorX();
0167 clustersOnDet.push_back(cluster);
0168 }
0169 }
0170 }
0171
0172
0173 theSeeds.clear();
0174
0175 for (auto itr = range.first; itr != range.second; ++itr) {
0176 const unsigned hitidx = itr->second;
0177 clear_buffer(begin + hitidx);
0178 }
0179 }
0180 }
0181
0182
0183
0184
0185
0186
0187
0188
0189 void MTDThresholdClusterizer::clear_buffer(RecHitIterator itr) { theBuffer.clear(itr->row(), itr->column()); }
0190
0191
0192
0193
0194 void MTDThresholdClusterizer::copy_to_buffer(RecHitIterator itr, const MTDGeometry* geom, const MTDTopology* topo) {
0195 MTDDetId mtdId = MTDDetId(itr->detid());
0196 int row = itr->row();
0197 int col = itr->column();
0198 GeomDetEnumerators::Location subDet = GeomDetEnumerators::invalidLoc;
0199 float energy = itr->energy();
0200 float time = itr->time();
0201 float timeError = itr->timeError();
0202 float position = itr->position();
0203 float xpos = 0.f;
0204
0205 LocalPoint local_point(0, 0, 0);
0206 LocalError local_error(0, 0, 0);
0207 if (mtdId.mtdSubDetector() == MTDDetId::BTL) {
0208 subDet = GeomDetEnumerators::barrel;
0209 BTLDetId id = itr->id();
0210 DetId geoId = id.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topo->getMTDTopologyMode()));
0211 const auto& det = geom->idToDet(geoId);
0212 const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(det->topology());
0213 const RectangularMTDTopology& topol = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0214
0215 LocalPoint lp_pixel(position, 0, 0);
0216 local_point = topol.pixelToModuleLocalPoint(lp_pixel, row, col);
0217 BTLRecHitsErrorEstimatorIM btlError(det, local_point);
0218 local_error = btlError.localError();
0219 xpos = local_point.x();
0220 } else if (mtdId.mtdSubDetector() == MTDDetId::ETL) {
0221 subDet = GeomDetEnumerators::endcap;
0222 ETLDetId id = itr->id();
0223 DetId geoId = id.geographicalId();
0224 const auto& det = geom->idToDet(geoId);
0225 const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(det->topology());
0226 const RectangularMTDTopology& topol = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0227
0228 LocalPoint lp_pixel(0, 0, 0);
0229 local_point = topol.pixelToModuleLocalPoint(lp_pixel, row, col);
0230 }
0231
0232 LogDebug("MTDThresholdClusterizer") << "DetId " << mtdId.rawId() << " subd " << mtdId.mtdSubDetector() << " row/col "
0233 << row << " / " << col << " energy " << energy << " time " << time
0234 << " time error " << timeError << " local_point " << local_point.x() << " "
0235 << local_point.y() << " " << local_point.z() << " local error "
0236 << std::sqrt(local_error.xx()) << " " << std::sqrt(local_error.yy()) << " xpos "
0237 << xpos;
0238 if (energy > theHitThreshold) {
0239 theBuffer.set(row, col, subDet, energy, time, timeError, local_error, local_point, xpos);
0240 if (energy > theSeedThreshold)
0241 theSeeds.push_back(FTLCluster::FTLHitPos(row, col));
0242
0243 }
0244 }
0245
0246
0247
0248
0249 FTLCluster MTDThresholdClusterizer::make_cluster(const FTLCluster::FTLHitPos& hit) {
0250
0251
0252 GeomDetEnumerators::Location seed_subdet = theBuffer.subDet(hit.row(), hit.col());
0253 float seed_energy = theBuffer.energy(hit.row(), hit.col());
0254 float seed_time = theBuffer.time(hit.row(), hit.col());
0255 float seed_time_error = theBuffer.time_error(hit.row(), hit.col());
0256 auto const seedPoint = theBuffer.local_point(hit.row(), hit.col());
0257 double seed_error_xx = theBuffer.local_error(hit.row(), hit.col()).xx();
0258 double seed_error_yy = theBuffer.local_error(hit.row(), hit.col()).yy();
0259 double seed_xpos = theBuffer.xpos(hit.row(), hit.col());
0260 theBuffer.clear(hit);
0261
0262 AccretionCluster acluster;
0263 acluster.add(hit, seed_energy, seed_time, seed_time_error);
0264
0265
0266 std::array<float, AccretionCluster::MAXSIZE> pixel_x{{-99999.}};
0267 std::array<float, AccretionCluster::MAXSIZE> pixel_errx2{{-99999.}};
0268 if (seed_subdet == GeomDetEnumerators::barrel) {
0269 pixel_x[acluster.top()] = seed_xpos;
0270 pixel_errx2[acluster.top()] = seed_error_xx;
0271 }
0272
0273 bool stopClus = false;
0274
0275 while (!acluster.empty() && !stopClus) {
0276
0277 auto curInd = acluster.top();
0278 acluster.pop();
0279 for (auto c = std::max(0, int(acluster.y[curInd] - 1));
0280 c < std::min(int(acluster.y[curInd] + 2), int(theBuffer.columns())) && !stopClus;
0281 ++c) {
0282 for (auto r = std::max(0, int(acluster.x[curInd] - 1));
0283 r < std::min(int(acluster.x[curInd] + 2), int(theBuffer.rows())) && !stopClus;
0284 ++r) {
0285 LogDebug("MTDThresholdClusterizer")
0286 << "Clustering subdet " << seed_subdet << " around " << curInd << " X,Y = " << acluster.x[curInd] << ","
0287 << acluster.y[curInd] << " r,c = " << r << "," << c << " energy,time = " << theBuffer.energy(r, c) << " "
0288 << theBuffer.time(r, c);
0289 if (theBuffer.energy(r, c) > theHitThreshold) {
0290 if (std::abs(theBuffer.time(r, c) - seed_time) >
0291 theTimeThreshold *
0292 sqrt(theBuffer.time_error(r, c) * theBuffer.time_error(r, c) + seed_time_error * seed_time_error))
0293 continue;
0294 if ((seed_subdet == GeomDetEnumerators::barrel) && (theBuffer.subDet(r, c) == GeomDetEnumerators::barrel) &&
0295 (thePositionThreshold > 0)) {
0296 double hit_error_xx = theBuffer.local_error(r, c).xx();
0297 double hit_error_yy = theBuffer.local_error(r, c).yy();
0298 if (((theBuffer.local_point(r, c) - seedPoint).mag2()) >
0299 thePositionThreshold * thePositionThreshold *
0300 (hit_error_xx + seed_error_xx + hit_error_yy + seed_error_yy))
0301 continue;
0302 }
0303 FTLCluster::FTLHitPos newhit(r, c);
0304 if (!acluster.add(newhit, theBuffer.energy(r, c), theBuffer.time(r, c), theBuffer.time_error(r, c))) {
0305 stopClus = true;
0306 break;
0307 }
0308 if (theBuffer.subDet(r, c) == GeomDetEnumerators::barrel) {
0309 pixel_x[acluster.top()] = theBuffer.xpos(r, c);
0310 pixel_errx2[acluster.top()] = theBuffer.local_error(r, c).xx();
0311 }
0312 theBuffer.clear(newhit);
0313 }
0314 }
0315 }
0316 }
0317
0318 FTLCluster cluster(theCurrentId,
0319 acluster.isize,
0320 acluster.energy.data(),
0321 acluster.time.data(),
0322 acluster.timeError.data(),
0323 acluster.x.data(),
0324 acluster.y.data(),
0325 acluster.xmin,
0326 acluster.ymin);
0327
0328
0329
0330 if (seed_subdet == GeomDetEnumerators::barrel) {
0331 float sumW(0.f), sumXW(0.f), sumXW2(0.f);
0332 for (unsigned int index = 0; index < acluster.top(); index++) {
0333 sumW += acluster.energy[index];
0334 sumXW += acluster.energy[index] * pixel_x[index];
0335 sumXW2 += acluster.energy[index] * acluster.energy[index] * pixel_errx2[index];
0336 }
0337 cluster.setClusterPosX(sumXW / sumW);
0338 cluster.setClusterErrorX(std::sqrt(sumXW2) / sumW);
0339 }
0340
0341 return cluster;
0342 }