DetLayers/src/MTDDetSector.cc

0001 //#define EDM_ML_DEBUG
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   // Add here the sector build based on a collection of GeomDets, mimic what done in ForwardDetRingOneZ
0035   // using the code from tracker BladeShapeBuilderFromDet
0036   // simple initial version, no sorting for the time being
0037   setDisk(MTDDiskSectorBuilderFromDet()(theDets));
0038 }
0039
0040 const vector<const GeometricSearchDet*>& MTDDetSector::components() const {
0041   // FIXME dummy impl.
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   // Propagate and check that the result is within bounds
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
0091   LogTrace("MTDDetLayers") << "Starting position: " << startPos << " starting p/pT: " << tsos.globalMomentum().mag()
0092                            << " / " << tsos.globalMomentum().perp();
0093
0094   // determine distance of det center from extrapolation on the surface, sort dets accordingly
0095
0096   size_t idetMin = basicComponents().size();
0097   double dist2Min = std::numeric_limits<double>::max();
0098   std::vector<std::pair<double, size_t> > tmpDets;
0099   tmpDets.reserve(basicComponents().size());
0100
0101   for (size_t idet = 0; idet < basicComponents().size(); idet++) {
0102     double dist2 = (startPos - theDets[idet]->position()).mag2();
0103     tmpDets.emplace_back(dist2, idet);
0104     if (dist2 < dist2Min) {
0105       dist2Min = dist2;
0106       idetMin = idet;
0107     }
0108   }
0109
0110   //look for the compatibledets considering each line of the sector
0111
0112   if (add(idetMin, result, tsos, prop, est)) {
0113     compatibleDetsLine(idetMin, result, tsos, prop, est);
0114
0115     for (int iside = -1; iside <= 1; iside += 2) {
0116       bool isCompatible(true);
0117       size_t idetNew(idetMin);
0118       size_t closest = theDets.size();
0119
0120       while (isCompatible) {
0121         idetNew = vshift(theDets[idetNew]->geographicalId().rawId(), iside, closest);
0122         if (idetNew >= theDets.size()) {
0123           if (closest < theDets.size()) {
0124             idetNew = closest;
0125           } else {
0126             break;
0127           }
0128         }
0129         isCompatible = add(idetNew, result, tsos, prop, est);
0130         if (isCompatible) {
0131           compatibleDetsLine(idetNew, result, tsos, prop, est);
0132         }
0133       }
0134     }
0135   }
0136
0137 #ifdef EDM_ML_DEBUG
0138   if (result.empty()) {
0139     LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets, closest not compatible!";
0140   } else {
0141     LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets, found " << result.size() << " compatible dets";
0142   }
0143 #endif
0144
0145   return result;
0146 }
0147
0148 void MTDDetSector::compatibleDetsV(const TrajectoryStateOnSurface&,
0149                                    const Propagator&,
0150                                    const MeasurementEstimator&,
0151                                    std::vector<DetWithState>&) const {
0152   edm::LogError("MTDDetLayers") << "At the moment not a real implementation";
0153 }
0154
0155 vector<DetGroup> MTDDetSector::groupedCompatibleDets(const TrajectoryStateOnSurface& startingState,
0156                                                      const Propagator& prop,
0157                                                      const MeasurementEstimator& est) const {
0158   // FIXME should be implemented to allow returning  overlapping chambers
0159   // as separate groups!
0160   edm::LogInfo("MTDDetLayers") << "dummy implementation of MTDDetSector::groupedCompatibleDets()";
0161   vector<DetGroup> result;
0162   return result;
0163 }
0164
0165 bool MTDDetSector::add(size_t idet,
0166                        vector<DetWithState>& result,
0167                        const TrajectoryStateOnSurface& tsos,
0168                        const Propagator& prop,
0169                        const MeasurementEstimator& est) const {
0170   pair<bool, TrajectoryStateOnSurface> compat = theCompatibilityChecker.isCompatible(theDets[idet], tsos, prop, est);
0171
0172   if (compat.first) {
0173     result.push_back(DetWithState(theDets[idet], compat.second));
0174     LogTrace("MTDDetLayers") << "MTDDetSector::compatibleDets found compatible det idetMin " << idet
0175                              << " detId = " << theDets[idet]->geographicalId().rawId() << " at "
0176                              << theDets[idet]->position()
0177                              << " dist = " << std::sqrt((tsos.globalPosition() - theDets[idet]->position()).mag2());
0178   }
0179
0180   return compat.first;
0181 }
0182
0183 std::ostream& operator<<(std::ostream& os, const MTDDetSector& id) {
0184   auto fround = [&](double in) {
0185     std::stringstream ss;
0186     ss << std::fixed << std::setprecision(3) << std::setw(14) << roundIfNear0(in);
0187     return ss.str();
0188   };
0189
0190   os << " MTDDetSector at " << std::fixed << std::setprecision(3) << roundVecIfNear0(id.specificSurface().position())
0191      << std::endl
0192      << " L/W/T   : " << fround(id.specificSurface().bounds().length()) << " / "
0193      << fround(id.specificSurface().bounds().width()) << " / " << fround(id.specificSurface().bounds().thickness())
0194      << std::endl
0195      << " rmin    : " << fround(id.specificSurface().innerRadius()) << std::endl
0196      << " rmax    : " << fround(id.specificSurface().outerRadius()) << std::endl
0197      << " phi ref : " << fround(id.specificSurface().position().phi()) << std::endl
0198      << " phi w/2 : " << fround(id.specificSurface().phiHalfExtension()) << std::endl;
0199   return os;
0200 }
0201
0202 void MTDDetSector::compatibleDetsLine(const size_t idetMin,
0203                                       vector<DetWithState>& result,
0204                                       const TrajectoryStateOnSurface& tsos,
0205                                       const Propagator& prop,
0206                                       const MeasurementEstimator& est) const {
0207   for (int iside = -1; iside <= 1; iside += 2) {
0208     bool isCompatible(true);
0209     size_t idetNew(idetMin);
0210
0211     while (isCompatible) {
0212       idetNew = hshift(theDets[idetNew]->geographicalId().rawId(), iside);
0213       if (idetNew >= theDets.size()) {
0214         break;
0215       }
0216       isCompatible = add(idetNew, result, tsos, prop, est);
0217     }
0218   }
0219
0220   return;
0221 }
0222
0223 size_t MTDDetSector::hshift(const uint32_t detid, const int horizontalShift) const {
0224   return topo_->hshiftETL(detid, horizontalShift);
0225 }
0226
0227 size_t MTDDetSector::vshift(const uint32_t detid, const int verticalShift, size_t& closest) const {
0228   return topo_->vshiftETL(detid, verticalShift, closest);
0229 }