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   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   // determine distance of det center from extrapolation on the surface, sort dets accordingly
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   //look for the compatibledets considering each line of the sector
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   // FIXME should be implemented to allow returning  overlapping chambers
0166   // as separate groups!
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 }