File indexing completed on 2022-05-26 22:39:47
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 + 1, ncols + 1);
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:" << cluster.time() << " X:" << cluster.x() << " Y:" << cluster.y();
0165 clustersOnDet.push_back(cluster);
0166 }
0167 }
0168 }
0169
0170
0171 theSeeds.clear();
0172
0173 for (auto itr = range.first; itr != range.second; ++itr) {
0174 const unsigned hitidx = itr->second;
0175 clear_buffer(begin + hitidx);
0176 }
0177 }
0178 }
0179
0180
0181
0182
0183
0184
0185
0186
0187 void MTDThresholdClusterizer::clear_buffer(RecHitIterator itr) { theBuffer.clear(itr->row(), itr->column()); }
0188
0189
0190
0191
0192 void MTDThresholdClusterizer::copy_to_buffer(RecHitIterator itr, const MTDGeometry* geom, const MTDTopology* topo) {
0193 MTDDetId mtdId = MTDDetId(itr->detid());
0194 int row = itr->row();
0195 int col = itr->column();
0196 GeomDetEnumerators::Location subDet = GeomDetEnumerators::invalidLoc;
0197 float energy = itr->energy();
0198 float time = itr->time();
0199 float timeError = itr->timeError();
0200 float position = itr->position();
0201
0202 LocalError local_error(0, 0, 0);
0203 GlobalPoint global_point(0, 0, 0);
0204 if (mtdId.mtdSubDetector() == MTDDetId::BTL) {
0205 subDet = GeomDetEnumerators::barrel;
0206 BTLDetId id = itr->id();
0207 DetId geoId = id.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topo->getMTDTopologyMode()));
0208 const auto& det = geom->idToDet(geoId);
0209 const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(det->topology());
0210 const RectangularMTDTopology& topol = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0211 MeasurementPoint mp(row, col);
0212 LocalPoint lp_ctr = topol.localPosition(mp);
0213 LocalPoint lp(lp_ctr.x() + position + topol.pitch().first * 0.5f, lp_ctr.y(), lp_ctr.z());
0214
0215
0216 global_point = det->toGlobal(lp);
0217 BTLRecHitsErrorEstimatorIM btlError(det, lp);
0218 local_error = btlError.localError();
0219 } else if (mtdId.mtdSubDetector() == MTDDetId::ETL) {
0220 subDet = GeomDetEnumerators::endcap;
0221 ETLDetId id = itr->id();
0222 DetId geoId = id.geographicalId();
0223 const auto& det = geom->idToDet(geoId);
0224 const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(det->topology());
0225 const RectangularMTDTopology& topol = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0226
0227 MeasurementPoint mp(row, col);
0228 LocalPoint lp = topol.localPosition(mp);
0229 global_point = det->toGlobal(lp);
0230 }
0231
0232 LogDebug("MTDThresholdClusterizer") << "ROW " << row << " COL " << col << " ENERGY " << energy << " TIME " << time;
0233 if (energy > theHitThreshold) {
0234 theBuffer.set(row, col, subDet, energy, time, timeError, local_error, global_point);
0235 if (energy > theSeedThreshold)
0236 theSeeds.push_back(FTLCluster::FTLHitPos(row, col));
0237
0238 }
0239 }
0240
0241
0242
0243
0244 FTLCluster MTDThresholdClusterizer::make_cluster(const FTLCluster::FTLHitPos& hit) {
0245
0246
0247 GeomDetEnumerators::Location seed_subdet = theBuffer.subDet(hit.row(), hit.col());
0248 float seed_energy = theBuffer.energy(hit.row(), hit.col());
0249 float seed_time = theBuffer.time(hit.row(), hit.col());
0250 float seed_time_error = theBuffer.time_error(hit.row(), hit.col());
0251 auto const seedPoint = theBuffer.global_point(hit.row(), hit.col());
0252 double seed_error_xx = theBuffer.local_error(hit.row(), hit.col()).xx();
0253 double seed_error_yy = theBuffer.local_error(hit.row(), hit.col()).yy();
0254 theBuffer.clear(hit);
0255
0256 AccretionCluster acluster;
0257 acluster.add(hit, seed_energy, seed_time, seed_time_error);
0258
0259 bool stopClus = false;
0260
0261 while (!acluster.empty() && !stopClus) {
0262
0263 auto curInd = acluster.top();
0264 acluster.pop();
0265 for (auto c = std::max(0, int(acluster.y[curInd] - 1));
0266 c < std::min(int(acluster.y[curInd] + 2), int(theBuffer.columns())) && !stopClus;
0267 ++c) {
0268 for (auto r = std::max(0, int(acluster.x[curInd] - 1));
0269 r < std::min(int(acluster.x[curInd] + 2), int(theBuffer.rows())) && !stopClus;
0270 ++r) {
0271 if (theBuffer.energy(r, c) > theHitThreshold) {
0272 if (std::abs(theBuffer.time(r, c) - seed_time) >
0273 theTimeThreshold *
0274 sqrt(theBuffer.time_error(r, c) * theBuffer.time_error(r, c) + seed_time_error * seed_time_error))
0275 continue;
0276 if ((seed_subdet == GeomDetEnumerators::barrel) && (theBuffer.subDet(r, c) == GeomDetEnumerators::barrel)) {
0277 double hit_error_xx = theBuffer.local_error(r, c).xx();
0278 double hit_error_yy = theBuffer.local_error(r, c).yy();
0279 if (thePositionThreshold > 0) {
0280 if (((theBuffer.global_point(r, c) - seedPoint).mag2()) >
0281 thePositionThreshold * thePositionThreshold *
0282 (hit_error_xx + seed_error_xx + hit_error_yy + seed_error_yy))
0283 continue;
0284 }
0285 }
0286 FTLCluster::FTLHitPos newhit(r, c);
0287 if (!acluster.add(newhit, theBuffer.energy(r, c), theBuffer.time(r, c), theBuffer.time_error(r, c))) {
0288 stopClus = true;
0289 break;
0290 }
0291 theBuffer.clear(newhit);
0292 }
0293 }
0294 }
0295 }
0296
0297 FTLCluster cluster(theCurrentId,
0298 acluster.isize,
0299 acluster.energy.data(),
0300 acluster.time.data(),
0301 acluster.timeError.data(),
0302 acluster.x.data(),
0303 acluster.y.data(),
0304 acluster.xmin,
0305 acluster.ymin);
0306 return cluster;
0307 }