DetLayers/src/MTDDetTray.cc

0001 //#define EDM_ML_DEBUG
0002
0003 /** \file
0004  *
0005  *  \author L. Gray - FNAL
0006  */
0007
0008 #include "RecoMTD/DetLayers/interface/MTDDetTray.h"
0009 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0010 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0011 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 #include <iostream>
0015
0016 #include "DataFormats/Math/interface/Rounding.h"
0017
0018 using namespace std;
0019 using namespace cms_rounding;
0020
0021 MTDDetTray::MTDDetTray(vector<const GeomDet*>::const_iterator first, vector<const GeomDet*>::const_iterator last)
0022     : DetRodOneR(first, last) {
0023   init();
0024 }
0025
0026 MTDDetTray::MTDDetTray(const vector<const GeomDet*>& vdets) : DetRodOneR(vdets) { init(); }
0027
0028 void MTDDetTray::init() { theBinFinder = BinFinderType(basicComponents().begin(), basicComponents().end()); }
0029
0030 MTDDetTray::~MTDDetTray() {}
0031
0032 const vector<const GeometricSearchDet*>& MTDDetTray::components() const {
0033   // FIXME dummy impl.
0034   edm::LogError("MTDDetLayers") << "temporary dummy implementation of MTDDetTray::components()!!";
0035   static const vector<const GeometricSearchDet*> result;
0036   return result;
0037 }
0038
0039 pair<bool, TrajectoryStateOnSurface> MTDDetTray::compatible(const TrajectoryStateOnSurface& ts,
0040                                                             const Propagator& prop,
0041                                                             const MeasurementEstimator& est) const {
0042   TrajectoryStateOnSurface ms = prop.propagate(ts, specificSurface());
0043   return make_pair(ms.isValid() and est.estimate(ms, specificSurface()) != 0, ms);
0044 }
0045
0046 vector<GeometricSearchDet::DetWithState> MTDDetTray::compatibleDets(const TrajectoryStateOnSurface& startingState,
0047                                                                     const Propagator& prop,
0048                                                                     const MeasurementEstimator& est) const {
0049   LogTrace("MTDDetLayers") << "MTDDetTray::compatibleDets, Surface at R,phi: " << surface().position().perp() << ","
0050                            << surface().position().phi() << "     DetRod pos."
0051                            << " TS at R,phi: " << startingState.globalPosition().perp() << ","
0052                            << startingState.globalPosition().phi() << " TSOS dx/dy "
0053                            << std::sqrt(startingState.localError().positionError().xx()) << " "
0054                            << std::sqrt(startingState.localError().positionError().yy());
0055
0056   vector<DetWithState> result;
0057
0058   // Propagate and check that the result is within bounds
0059   pair<bool, TrajectoryStateOnSurface> compat = compatible(startingState, prop, est);
0060
0061   if (!compat.first) {
0062     LogTrace("MTDDetLayers") << "    MTDDetTray::compatibleDets: not compatible"
0063                              << "    (should not have been selected!)";
0064     return result;
0065   }
0066
0067   // Find the most probable destination component
0068   TrajectoryStateOnSurface& tsos = compat.second;
0069   GlobalPoint startPos = tsos.globalPosition();
0070   int closest = theBinFinder.binIndex(startPos.z());
0071   const vector<const GeomDet*> dets = basicComponents();
0072   LogTrace("MTDDetLayers") << "     MTDDetTray::compatibleDets, closest det: " << closest
0073                            << " pos: " << dets[closest]->surface().position() << " impact " << startPos;
0074
0075   // Add this detector, if it is compatible
0076   // NOTE: add performs a null propagation
0077   add(closest, result, tsos, prop, est);
0078
0079 #ifdef EDM_ML_DEBUG
0080   int nclosest = result.size();
0081   int nnextdet = 0;  // just DEBUG counters
0082 #endif
0083
0084   // Try the neighbors on each side until no more compatible.
0085   // If closest is not compatible the next cannot be either
0086   if (!result.empty()) {
0087     const BoundPlane& closestPlane(dets[closest]->surface());
0088     MeasurementEstimator::Local2DVector maxDistance = est.maximalLocalDisplacement(result.front().second, closestPlane);
0089
0090     // detHalfLen is assumed to be the same for all detectors.
0091     float detHalfLen = closestPlane.bounds().length() / 2.;
0092
0093     for (unsigned int idet = closest + 1; idet < dets.size(); idet++) {
0094       LocalPoint nextPos(dets[idet]->toLocal(startPos));
0095       if (fabs(nextPos.y()) < detHalfLen + maxDistance.y()) {
0096 #ifdef EDM_ML_DEBUG
0097         LogTrace("MTDDetLayers") << "     negativeZ: det:" << idet << " pos " << nextPos.y() << " maxDistance "
0098                                  << maxDistance.y();
0099         nnextdet++;
0100 #endif
0101         if (!add(idet, result, tsos, prop, est))
0102           break;
0103       } else {
0104         break;
0105       }
0106     }
0107
0108     for (int idet = closest - 1; idet >= 0; idet--) {
0109       LocalPoint nextPos(dets[idet]->toLocal(startPos));
0110       if (fabs(nextPos.y()) < detHalfLen + maxDistance.y()) {
0111 #ifdef EDM_ML_DEBUG
0112         LogTrace("MTDDetLayers") << "     positiveZ: det:" << idet << " pos " << nextPos.y() << " maxDistance "
0113                                  << maxDistance.y();
0114         nnextdet++;
0115 #endif
0116         if (!add(idet, result, tsos, prop, est))
0117           break;
0118       } else {
0119         break;
0120       }
0121     }
0122   }
0123
0124 #ifdef EDM_ML_DEBUG
0125   LogTrace("MTDDetLayers") << "     MTDDetTray::compatibleDets, size: " << result.size() << " on closest: " << nclosest
0126                            << " # checked dets: " << nnextdet + 1;
0127   if (result.empty()) {
0128     LogTrace("MTDDetLayers") << "   ***Rod not compatible---should have been discarded before!!!";
0129   }
0130 #endif
0131   return result;
0132 }
0133
0134 vector<DetGroup> MTDDetTray::groupedCompatibleDets(const TrajectoryStateOnSurface& startingState,
0135                                                    const Propagator& prop,
0136                                                    const MeasurementEstimator& est) const {
0137   vector<GeometricSearchDet::DetWithState> detWithStates;
0138
0139   detWithStates = compatibleDets(startingState, prop, est);
0140
0141   vector<DetGroup> result;
0142   if (!detWithStates.empty()) {
0143     result.push_back(DetGroup(detWithStates));
0144   }
0145   LogTrace("MTDDetLayers") << "MTDdetTray Compatible modules: " << result.size();
0146   return result;
0147 }