File indexing completed on 2025-04-06 22:43:12
0001
0002
0003 #include "RecoMTD/DetLayers/interface/MTDDetSector.h"
0004 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0005 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0006 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0007 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0008 #include "MTDDiskSectorBuilderFromDet.h"
0009 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 #include <iostream>
0013 #include <iomanip>
0014 #include <vector>
0015
0016 #include "DataFormats/Math/interface/Rounding.h"
0017
0018 using namespace std;
0019 using namespace cms_rounding;
0020
0021 MTDDetSector::MTDDetSector(vector<const GeomDet*>::const_iterator first,
0022 vector<const GeomDet*>::const_iterator last,
0023 const MTDTopology& topo)
0024 : GeometricSearchDet(false), theDets(first, last), topo_(&topo) {
0025 init();
0026 }
0027
0028 MTDDetSector::MTDDetSector(const vector<const GeomDet*>& vdets, const MTDTopology& topo)
0029 : GeometricSearchDet(false), theDets(vdets), topo_(&topo) {
0030 init();
0031 }
0032
0033 void MTDDetSector::init() {
0034
0035
0036
0037 setDisk(MTDDiskSectorBuilderFromDet()(theDets));
0038 }
0039
0040 const vector<const GeometricSearchDet*>& MTDDetSector::components() const {
0041
0042 edm::LogError("MTDDetLayers") << "temporary dummy implementation of MTDDetSector::components()!!";
0043 static const vector<const GeometricSearchDet*> result;
0044 return result;
0045 }
0046
0047 pair<bool, TrajectoryStateOnSurface> MTDDetSector::compatible(const TrajectoryStateOnSurface& ts,
0048 const Propagator& prop,
0049 const MeasurementEstimator& est) const {
0050 TrajectoryStateOnSurface ms = prop.propagate(ts, specificSurface());
0051
0052 #ifdef EDM_ML_DEBUG
0053 LogTrace("MTDDetLayers") << "MTDDetSector::compatible, sector: \n"
0054 << (*this) << "\n TS at Z,R,phi: " << std::fixed << std::setw(14) << ts.globalPosition().z()
0055 << " , " << std::setw(14) << ts.globalPosition().perp() << " , " << std::setw(14)
0056 << ts.globalPosition().phi();
0057 if (ms.isValid()) {
0058 LogTrace("MTDDetLayers") << " DEST at Z,R,phi: " << std::fixed << std::setw(14) << ms.globalPosition().z() << " , "
0059 << std::setw(14) << ms.globalPosition().perp() << " , " << std::setw(14)
0060 << ms.globalPosition().phi() << " local Z: " << std::setw(14) << ms.localPosition().z();
0061 } else {
0062 LogTrace("MTDDetLayers") << " DEST: not valid";
0063 }
0064 #endif
0065
0066 return make_pair(ms.isValid() and est.estimate(ms, specificSurface()) != 0, ms);
0067 }
0068
0069 vector<GeometricSearchDet::DetWithState> MTDDetSector::compatibleDets(const TrajectoryStateOnSurface& startingState,
0070 const Propagator& prop,
0071 const MeasurementEstimator& est) const {
0072 LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets, sector: \n"
0073 << (*this) << "\n TS at Z,R,phi: " << std::fixed << std::setw(14)
0074 << startingState.globalPosition().z() << " , " << std::setw(14)
0075 << startingState.globalPosition().perp() << " , " << std::setw(14)
0076 << startingState.globalPosition().phi();
0077
0078 vector<DetWithState> result;
0079
0080
0081 pair<bool, TrajectoryStateOnSurface> compat = compatible(startingState, prop, est);
0082 if (!compat.first) {
0083 LogTrace("MTDDetLayers") << " MTDDetSector::compatibleDets: not compatible"
0084 << " (should not have been selected!)";
0085 return result;
0086 }
0087
0088 TrajectoryStateOnSurface& tsos = compat.second;
0089 GlobalPoint startPos = tsos.globalPosition();
0090 LocalTrajectoryError startLocalErr = tsos.localError();
0091
0092 LogTrace("MTDDetLayers") << "Starting position: " << startPos << " starting p/pT: " << tsos.globalMomentum().mag()
0093 << " / " << tsos.globalMomentum().perp() << " local error x/y "
0094 << std::sqrt(startLocalErr.positionError().xx()) << " "
0095 << std::sqrt(startLocalErr.positionError().yy());
0096
0097
0098
0099 size_t idetMin = basicComponents().size();
0100 double dist2Min = std::numeric_limits<double>::max();
0101 std::vector<std::pair<double, size_t> > tmpDets;
0102 tmpDets.reserve(basicComponents().size());
0103
0104 for (size_t idet = 0; idet < basicComponents().size(); idet++) {
0105 double dist2 = (startPos - theDets[idet]->position()).mag2();
0106 tmpDets.emplace_back(dist2, idet);
0107 if (dist2 < dist2Min) {
0108 dist2Min = dist2;
0109 idetMin = idet;
0110 }
0111 }
0112
0113
0114
0115 if (add(idetMin, result, tsos, prop, est)) {
0116 compatibleDetsLine(idetMin, result, tsos, prop, est);
0117
0118 for (int iside = -1; iside <= 1; iside += 2) {
0119 bool isCompatible(true);
0120 size_t idetNew(idetMin);
0121 size_t closest = theDets.size();
0122 LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets, maximum size available " << closest;
0123
0124 while (isCompatible) {
0125 idetNew = vshift(theDets[idetNew]->geographicalId().rawId(), iside, closest);
0126 LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets, vshift iside " << iside << " idetNew " << idetNew
0127 << " closest " << closest;
0128 if (idetNew >= theDets.size()) {
0129 if (closest < theDets.size()) {
0130 idetNew = closest;
0131 } else {
0132 break;
0133 }
0134 }
0135 isCompatible = add(idetNew, result, tsos, prop, est);
0136 if (isCompatible) {
0137 compatibleDetsLine(idetNew, result, tsos, prop, est);
0138 }
0139 }
0140 LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets, exiting loop on vshift iside " << iside;
0141 }
0142 }
0143
0144 #ifdef EDM_ML_DEBUG
0145 if (result.empty()) {
0146 LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets, closest not compatible!";
0147 } else {
0148 LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets, found " << result.size() << " compatible dets";
0149 }
0150 #endif
0151
0152 return result;
0153 }
0154
0155 void MTDDetSector::compatibleDetsV(const TrajectoryStateOnSurface&,
0156 const Propagator&,
0157 const MeasurementEstimator&,
0158 std::vector<DetWithState>&) const {
0159 edm::LogError("MTDDetLayers") << "At the moment not a real implementation";
0160 }
0161
0162 vector<DetGroup> MTDDetSector::groupedCompatibleDets(const TrajectoryStateOnSurface& startingState,
0163 const Propagator& prop,
0164 const MeasurementEstimator& est) const {
0165
0166
0167 edm::LogInfo("MTDDetLayers") << "dummy implementation of MTDDetSector::groupedCompatibleDets()";
0168 vector<DetGroup> result;
0169 return result;
0170 }
0171
0172 bool MTDDetSector::add(size_t idet,
0173 vector<DetWithState>& result,
0174 const TrajectoryStateOnSurface& tsos,
0175 const Propagator& prop,
0176 const MeasurementEstimator& est) const {
0177 pair<bool, TrajectoryStateOnSurface> compat = theCompatibilityChecker.isCompatible(theDets[idet], tsos, prop, est);
0178
0179 if (compat.first) {
0180 result.push_back(DetWithState(theDets[idet], compat.second));
0181 LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets found compatible det idetMin " << idet + 1
0182 << " detId = " << theDets[idet]->geographicalId().rawId() << " at "
0183 << theDets[idet]->position()
0184 << " dist = " << std::sqrt((tsos.globalPosition() - theDets[idet]->position()).mag2());
0185 }
0186 if (result.size() > basicComponents().size()) {
0187 throw cms::Exception("MTDDetLayersFailure")
0188 << "ETL compatibleDets in sector in excess of collection size: " << basicComponents().size();
0189 }
0190
0191 return compat.first;
0192 }
0193
0194 std::ostream& operator<<(std::ostream& os, const MTDDetSector& id) {
0195 auto fround = [&](double in) {
0196 std::stringstream ss;
0197 ss << std::fixed << std::setprecision(3) << std::setw(14) << roundIfNear0(in);
0198 return ss.str();
0199 };
0200
0201 os << " MTDDetSector at " << std::fixed << std::setprecision(3) << roundVecIfNear0(id.specificSurface().position())
0202 << std::endl
0203 << " L/W/T : " << fround(id.specificSurface().bounds().length()) << " / "
0204 << fround(id.specificSurface().bounds().width()) << " / " << fround(id.specificSurface().bounds().thickness())
0205 << std::endl
0206 << " rmin : " << fround(id.specificSurface().innerRadius()) << std::endl
0207 << " rmax : " << fround(id.specificSurface().outerRadius()) << std::endl
0208 << " phi ref : " << fround(id.specificSurface().position().phi()) << std::endl
0209 << " phi w/2 : " << fround(id.specificSurface().phiHalfExtension()) << std::endl;
0210 return os;
0211 }
0212
0213 void MTDDetSector::compatibleDetsLine(const size_t idetMin,
0214 vector<DetWithState>& result,
0215 const TrajectoryStateOnSurface& tsos,
0216 const Propagator& prop,
0217 const MeasurementEstimator& est) const {
0218 for (int iside = -1; iside <= 1; iside += 2) {
0219 bool isCompatible(true);
0220 size_t idetNew(idetMin);
0221
0222 while (isCompatible) {
0223 idetNew = hshift(theDets[idetNew]->geographicalId().rawId(), iside);
0224 if (idetNew >= theDets.size()) {
0225 break;
0226 }
0227 isCompatible = add(idetNew, result, tsos, prop, est);
0228 }
0229 }
0230
0231 return;
0232 }
0233
0234 size_t MTDDetSector::hshift(const uint32_t detid, const int horizontalShift) const {
0235 return topo_->hshiftETL(detid, horizontalShift);
0236 }
0237
0238 size_t MTDDetSector::vshift(const uint32_t detid, const int verticalShift, size_t& closest) const {
0239 return topo_->vshiftETL(detid, verticalShift, closest);
0240 }