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