File indexing completed on 2023-03-17 11:20:16
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 using namespace std;
0017
0018 MTDDetSector::MTDDetSector(vector<const GeomDet*>::const_iterator first,
0019 vector<const GeomDet*>::const_iterator last,
0020 const MTDTopology& topo)
0021 : GeometricSearchDet(false), theDets(first, last), topo_(&topo) {
0022 init();
0023 }
0024
0025 MTDDetSector::MTDDetSector(const vector<const GeomDet*>& vdets, const MTDTopology& topo)
0026 : GeometricSearchDet(false), theDets(vdets), topo_(&topo) {
0027 init();
0028 }
0029
0030 void MTDDetSector::init() {
0031
0032
0033
0034 setDisk(MTDDiskSectorBuilderFromDet()(theDets));
0035 }
0036
0037 const vector<const GeometricSearchDet*>& MTDDetSector::components() const {
0038
0039 edm::LogError("MTDDetLayers") << "temporary dummy implementation of MTDDetSector::components()!!";
0040 static const vector<const GeometricSearchDet*> result;
0041 return result;
0042 }
0043
0044 pair<bool, TrajectoryStateOnSurface> MTDDetSector::compatible(const TrajectoryStateOnSurface& ts,
0045 const Propagator& prop,
0046 const MeasurementEstimator& est) const {
0047 TrajectoryStateOnSurface ms = prop.propagate(ts, specificSurface());
0048
0049 #ifdef EDM_ML_DEBUG
0050 LogTrace("MTDDetLayers") << "MTDDetSector::compatible, sector: \n"
0051 << (*this) << "\n TS at Z,R,phi: " << std::fixed << std::setw(14) << ts.globalPosition().z()
0052 << " , " << std::setw(14) << ts.globalPosition().perp() << " , " << std::setw(14)
0053 << ts.globalPosition().phi();
0054 if (ms.isValid()) {
0055 LogTrace("MTDDetLayers") << " DEST at Z,R,phi: " << std::fixed << std::setw(14) << ms.globalPosition().z() << " , "
0056 << std::setw(14) << ms.globalPosition().perp() << " , " << std::setw(14)
0057 << ms.globalPosition().phi() << " local Z: " << std::setw(14) << ms.localPosition().z();
0058 } else {
0059 LogTrace("MTDDetLayers") << " DEST: not valid";
0060 }
0061 #endif
0062
0063 return make_pair(ms.isValid() and est.estimate(ms, specificSurface()) != 0, ms);
0064 }
0065
0066 vector<GeometricSearchDet::DetWithState> MTDDetSector::compatibleDets(const TrajectoryStateOnSurface& startingState,
0067 const Propagator& prop,
0068 const MeasurementEstimator& est) const {
0069 LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets, sector: \n"
0070 << (*this) << "\n TS at Z,R,phi: " << std::fixed << std::setw(14)
0071 << startingState.globalPosition().z() << " , " << std::setw(14)
0072 << startingState.globalPosition().perp() << " , " << std::setw(14)
0073 << startingState.globalPosition().phi();
0074
0075 vector<DetWithState> result;
0076
0077
0078 pair<bool, TrajectoryStateOnSurface> compat = compatible(startingState, prop, est);
0079 if (!compat.first) {
0080 LogTrace("MTDDetLayers") << " MTDDetSector::compatibleDets: not compatible"
0081 << " (should not have been selected!)";
0082 return result;
0083 }
0084
0085 TrajectoryStateOnSurface& tsos = compat.second;
0086 GlobalPoint startPos = tsos.globalPosition();
0087
0088 LogTrace("MTDDetLayers") << "Starting position: " << startPos << " starting p/pT: " << tsos.globalMomentum().mag()
0089 << " / " << tsos.globalMomentum().perp();
0090
0091
0092
0093 size_t idetMin = basicComponents().size();
0094 double dist2Min = std::numeric_limits<double>::max();
0095 std::vector<std::pair<double, size_t> > tmpDets;
0096 tmpDets.reserve(basicComponents().size());
0097
0098 for (size_t idet = 0; idet < basicComponents().size(); idet++) {
0099 double dist2 = (startPos - theDets[idet]->position()).mag2();
0100 tmpDets.emplace_back(dist2, idet);
0101 if (dist2 < dist2Min) {
0102 dist2Min = dist2;
0103 idetMin = idet;
0104 }
0105 }
0106
0107
0108
0109 if (add(idetMin, result, tsos, prop, est)) {
0110 compatibleDetsLine(idetMin, result, tsos, prop, est);
0111
0112 for (int iside = -1; iside <= 1; iside += 2) {
0113 bool isCompatible(true);
0114 size_t idetNew(idetMin);
0115 size_t closest = theDets.size();
0116
0117 while (isCompatible) {
0118 idetNew = vshift(theDets[idetNew]->geographicalId().rawId(), iside, closest);
0119 if (idetNew >= theDets.size()) {
0120 if (closest < theDets.size()) {
0121 idetNew = closest;
0122 } else {
0123 break;
0124 }
0125 }
0126 isCompatible = add(idetNew, result, tsos, prop, est);
0127 if (isCompatible) {
0128 compatibleDetsLine(idetNew, result, tsos, prop, est);
0129 }
0130 }
0131 }
0132 }
0133
0134 #ifdef EDM_ML_DEBUG
0135 if (result.empty()) {
0136 LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets, closest not compatible!";
0137 } else {
0138 LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets, found " << result.size() << " compatible dets";
0139 }
0140 #endif
0141
0142 return result;
0143 }
0144
0145 void MTDDetSector::compatibleDetsV(const TrajectoryStateOnSurface&,
0146 const Propagator&,
0147 const MeasurementEstimator&,
0148 std::vector<DetWithState>&) const {
0149 edm::LogError("MTDDetLayers") << "At the moment not a real implementation";
0150 }
0151
0152 vector<DetGroup> MTDDetSector::groupedCompatibleDets(const TrajectoryStateOnSurface& startingState,
0153 const Propagator& prop,
0154 const MeasurementEstimator& est) const {
0155
0156
0157 edm::LogInfo("MTDDetLayers") << "dummy implementation of MTDDetSector::groupedCompatibleDets()";
0158 vector<DetGroup> result;
0159 return result;
0160 }
0161
0162 bool MTDDetSector::add(size_t idet,
0163 vector<DetWithState>& result,
0164 const TrajectoryStateOnSurface& tsos,
0165 const Propagator& prop,
0166 const MeasurementEstimator& est) const {
0167 pair<bool, TrajectoryStateOnSurface> compat = theCompatibilityChecker.isCompatible(theDets[idet], tsos, prop, est);
0168
0169 if (compat.first) {
0170 result.push_back(DetWithState(theDets[idet], compat.second));
0171 LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets found compatible det idetMin " << idet
0172 << " detId = " << theDets[idet]->geographicalId().rawId() << " at "
0173 << theDets[idet]->position()
0174 << " dist = " << std::sqrt((tsos.globalPosition() - theDets[idet]->position()).mag2());
0175 }
0176
0177 return compat.first;
0178 }
0179
0180 std::ostream& operator<<(std::ostream& os, const MTDDetSector& id) {
0181 os << " MTDDetSector at " << std::fixed << id.specificSurface().position() << std::endl
0182 << " L/W/T : " << std::setw(14) << id.specificSurface().bounds().length() << " / " << std::setw(14)
0183 << id.specificSurface().bounds().width() << " / " << std::setw(14) << id.specificSurface().bounds().thickness()
0184 << std::endl
0185 << " rmin : " << std::setw(14) << id.specificSurface().innerRadius() << std::endl
0186 << " rmax : " << std::setw(14) << id.specificSurface().outerRadius() << std::endl
0187 << " phi ref : " << std::setw(14) << id.specificSurface().position().phi() << std::endl
0188 << " phi w/2 : " << std::setw(14) << id.specificSurface().phiHalfExtension() << std::endl;
0189 return os;
0190 }
0191
0192 void MTDDetSector::compatibleDetsLine(const size_t idetMin,
0193 vector<DetWithState>& result,
0194 const TrajectoryStateOnSurface& tsos,
0195 const Propagator& prop,
0196 const MeasurementEstimator& est) const {
0197 for (int iside = -1; iside <= 1; iside += 2) {
0198 bool isCompatible(true);
0199 size_t idetNew(idetMin);
0200
0201 while (isCompatible) {
0202 idetNew = hshift(theDets[idetNew]->geographicalId().rawId(), iside);
0203 if (idetNew >= theDets.size()) {
0204 break;
0205 }
0206 isCompatible = add(idetNew, result, tsos, prop, est);
0207 }
0208 }
0209
0210 return;
0211 }
0212
0213 size_t MTDDetSector::hshift(const uint32_t detid, const int horizontalShift) const {
0214 return topo_->hshiftETL(detid, horizontalShift);
0215 }
0216
0217 size_t MTDDetSector::vshift(const uint32_t detid, const int verticalShift, size_t& closest) const {
0218 return topo_->vshiftETL(detid, verticalShift, closest);
0219 }