DetLayers/src/MTDRingForwardLayer.cc

0001 //#define EDM_ML_DEBUG
0002
0003 /** \file
0004  *
0005  *  \author L. Gray - FNAL
0006  */
0007
0008 #include <RecoMTD/DetLayers/interface/MTDRingForwardLayer.h>
0009 #include <RecoMTD/DetLayers/interface/MTDDetRing.h>
0010 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0011 #include <DataFormats/GeometrySurface/interface/SimpleDiskBounds.h>
0012 #include <TrackingTools/GeomPropagators/interface/Propagator.h>
0013 #include <TrackingTools/DetLayers/interface/MeasurementEstimator.h>
0014 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0015
0016 #include "TrackingTools/DetLayers/interface/RBorderFinder.h"
0017 #include "TrackingTools/DetLayers/interface/GeneralBinFinderInR.h"
0018
0019 #include <algorithm>
0020 #include <iostream>
0021 #include <vector>
0022
0023 using namespace std;
0024
0025 MTDRingForwardLayer::MTDRingForwardLayer(const vector<const ForwardDetRing*>& rings)
0026     : RingedForwardLayer(false),
0027       theRings(rings),
0028       theComponents(theRings.begin(), theRings.end()),
0029       theBinFinder(nullptr),
0030       isOverlapping(false) {
0031   // Initial values for R and Z bounds
0032   float theRmin = rings.front()->basicComponents().front()->position().perp();
0033   float theRmax = theRmin;
0034   float theZmin = rings.front()->position().z();
0035   float theZmax = theZmin;
0036
0037   // Cache chamber pointers (the basic components_)
0038   // and find extension in R and Z
0039   for (const auto& it : rings) {
0040     vector<const GeomDet*> tmp2 = it->basicComponents();
0041     theBasicComps.insert(theBasicComps.end(), tmp2.begin(), tmp2.end());
0042
0043     theRmin = min(theRmin, it->specificSurface().innerRadius());
0044     theRmax = max(theRmax, it->specificSurface().outerRadius());
0045     float halfThick = it->surface().bounds().thickness() / 2.;
0046     float zCenter = it->surface().position().z();
0047     theZmin = min(theZmin, zCenter - halfThick);
0048     theZmax = max(theZmax, zCenter + halfThick);
0049   }
0050
0051   RBorderFinder bf(theRings);
0052   isOverlapping = bf.isROverlapping();
0053   theBinFinder = new GeneralBinFinderInR<double>(bf);
0054
0055   // Build surface
0056
0057   float zPos = (theZmax + theZmin) / 2.;
0058   PositionType pos(0., 0., zPos);
0059   RotationType rot;
0060
0061   setSurface(new BoundDisk(pos, rot, new SimpleDiskBounds(theRmin, theRmax, theZmin - zPos, theZmax - zPos)));
0062
0063   LogTrace("MTDDetLayers") << "Constructing MTDRingForwardLayer: " << basicComponents().size() << " Dets "
0064                            << theRings.size() << " Rings "
0065                            << " Z: " << specificSurface().position().z() << " R1: " << specificSurface().innerRadius()
0066                            << " R2: " << specificSurface().outerRadius() << " Per.: " << bf.isRPeriodic()
0067                            << " Overl.: " << bf.isROverlapping();
0068 }
0069
0070 MTDRingForwardLayer::~MTDRingForwardLayer() {
0071   delete theBinFinder;
0072   for (auto& i : theRings) {
0073     delete i;
0074   }
0075 }
0076
0077 vector<GeometricSearchDet::DetWithState> MTDRingForwardLayer::compatibleDets(
0078     const TrajectoryStateOnSurface& startingState, const Propagator& prop, const MeasurementEstimator& est) const {
0079   vector<DetWithState> result;
0080
0081   LogTrace("MTDDetLayers") << "MTDRingForwardLayer::compatibleDets,"
0082                            << " R1 " << specificSurface().innerRadius() << " R2: " << specificSurface().outerRadius()
0083                            << " FTS at R: " << startingState.globalPosition().perp();
0084
0085   pair<bool, TrajectoryStateOnSurface> compat = compatible(startingState, prop, est);
0086
0087   if (!compat.first) {
0088     LogTrace("MTDDetLayers") << "     MTDRingForwardLayer::compatibleDets: not compatible"
0089                              << " (should not have been selected!)";
0090     return result;
0091   }
0092
0093   TrajectoryStateOnSurface& tsos = compat.second;
0094
0095   int closest = theBinFinder->binIndex(tsos.globalPosition().perp());
0096   const ForwardDetRing* closestRing = theRings[closest];
0097
0098   // Check the closest ring
0099
0100   LogTrace("MTDDetLayers") << "     MTDRingForwardLayer::fastCompatibleDets, closestRing: " << closest << " R1 "
0101                            << closestRing->specificSurface().innerRadius()
0102                            << " R2: " << closestRing->specificSurface().outerRadius()
0103                            << " FTS R: " << tsos.globalPosition().perp();
0104   if (tsos.hasError()) {
0105     LogTrace("MTDDetLayers") << " sR: " << sqrt(tsos.localError().positionError().yy())
0106                              << " sX: " << sqrt(tsos.localError().positionError().xx());
0107   }
0108
0109   result = closestRing->compatibleDets(tsos, prop, est);
0110
0111 #ifdef EDM_ML_DEBUG
0112   int nclosest = result.size();
0113   int nnextdet = 0;  // MDEBUG counters
0114 #endif
0115
0116   //FIXME: if closest is not compatible next cannot be either?
0117
0118   // Use state on layer surface. Note that local coordinates and errors
0119   // are the same on the layer and on all rings surfaces, since
0120   // all BoundDisks are centered in 0,0 and have the same rotation.
0121   // CAVEAT: if the rings are not at the same Z, the local position and error
0122   // will be "Z-projected" to the rings. This is a fairly good approximation.
0123   // However in this case additional propagation will be done when calling
0124   // compatibleDets.
0125   GlobalPoint startPos = tsos.globalPosition();
0126   LocalPoint nextPos(surface().toLocal(startPos));
0127
0128   for (unsigned int idet = closest + 1; idet < theRings.size(); idet++) {
0129     bool inside = false;
0130     if (tsos.hasError()) {
0131       inside = theRings[idet]->specificSurface().bounds().inside(nextPos, tsos.localError().positionError());
0132     } else {
0133       inside = theRings[idet]->specificSurface().bounds().inside(nextPos);
0134     }
0135     if (inside) {
0136 #ifdef EDM_ML_DEBUG
0137       LogTrace("MTDDetLayers") << "     MTDRingForwardLayer::fastCompatibleDets:NextRing" << idet << " R1 "
0138                                << theRings[idet]->specificSurface().innerRadius()
0139                                << " R2: " << theRings[idet]->specificSurface().outerRadius() << " FTS R "
0140                                << nextPos.perp();
0141       nnextdet++;
0142 #endif
0143       vector<DetWithState> nextRodDets = theRings[idet]->compatibleDets(tsos, prop, est);
0144       if (!nextRodDets.empty()) {
0145         result.insert(result.end(), nextRodDets.begin(), nextRodDets.end());
0146       } else {
0147         break;
0148       }
0149     }
0150   }
0151
0152   for (int idet = closest - 1; idet >= 0; idet--) {
0153     bool inside = false;
0154     if (tsos.hasError()) {
0155       inside = theRings[idet]->specificSurface().bounds().inside(nextPos, tsos.localError().positionError());
0156     } else {
0157       inside = theRings[idet]->specificSurface().bounds().inside(nextPos);
0158     }
0159     if (inside) {
0160 #ifdef EDM_ML_DEBUG
0161       LogTrace("MTDDetLayers") << "     MTDRingForwardLayer::fastCompatibleDets:PreviousRing:" << idet << " R1 "
0162                                << theRings[idet]->specificSurface().innerRadius()
0163                                << " R2: " << theRings[idet]->specificSurface().outerRadius() << " FTS R "
0164                                << nextPos.perp();
0165       nnextdet++;
0166 #endif
0167       vector<DetWithState> nextRodDets = theRings[idet]->compatibleDets(tsos, prop, est);
0168       if (!nextRodDets.empty()) {
0169         result.insert(result.end(), nextRodDets.begin(), nextRodDets.end());
0170       } else {
0171         break;
0172       }
0173     }
0174   }
0175
0176 #ifdef EDM_ML_DEBUG
0177   LogTrace("MTDDetLayers") << "     MTDRingForwardLayer::fastCompatibleDets: found: " << result.size()
0178                            << " on closest: " << nclosest << " # checked rings: " << 1 + nnextdet;
0179 #endif
0180
0181   return result;
0182 }
0183
0184 vector<DetGroup> MTDRingForwardLayer::groupedCompatibleDets(const TrajectoryStateOnSurface& startingState,
0185                                                             const Propagator& prop,
0186                                                             const MeasurementEstimator& est) const {
0187   // FIXME should return only 1 group
0188   edm::LogError("MTDDetLayers") << "dummy implementation of MTDRingForwardLayer::groupedCompatibleDets()";
0189   return vector<DetGroup>();
0190 }
0191
0192 GeomDetEnumerators::SubDetector MTDRingForwardLayer::subDetector() const {
0193   return theBasicComps.front()->subDetector();
0194 }
0195
0196 const vector<const GeometricSearchDet*>& MTDRingForwardLayer::components() const { return theComponents; }