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 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   // Add here the sector build based on a collection of GeomDets, mimic what done in ForwardDetRingOneZ
0032   // using the code from tracker BladeShapeBuilderFromDet
0033   // simple initial version, no sorting for the time being
0034   setDisk(MTDDiskSectorBuilderFromDet()(theDets));
0035 }
0036
0037 const vector<const GeometricSearchDet*>& MTDDetSector::components() const {
0038   // FIXME dummy impl.
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   // Propagate and check that the result is within bounds
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   // determine distance of det center from extrapolation on the surface, sort dets accordingly
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   //look for the compatibledets considering each line of the sector
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   // FIXME should be implemented to allow returning  overlapping chambers
0156   // as separate groups!
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 }