DetLayers/src/MTDTrayBarrelLayer.cc

0001 //#define EDM_ML_DEBUG
0002
0003 /** \file
0004  *
0005  *  \author L. Gray - FNAL
0006  */
0007
0008 #include "RecoMTD/DetLayers/interface/MTDTrayBarrelLayer.h"
0009 #include "RecoMTD/DetLayers/interface/MTDDetTray.h"
0010 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0011 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0012 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0013 #include "Utilities/BinningTools/interface/PeriodicBinFinderInPhi.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 #include <Utilities/General/interface/precomputed_value_sort.h>
0017 #include "DataFormats/GeometrySurface/interface/GeometricSorting.h"
0018
0019 #include "TrackingTools/DetLayers/interface/GeneralBinFinderInPhi.h"
0020 #include "TrackingTools/DetLayers/interface/PhiBorderFinder.h"
0021
0022 #include <algorithm>
0023 #include <iostream>
0024
0025 using namespace std;
0026
0027 MTDTrayBarrelLayer::MTDTrayBarrelLayer(vector<const DetRod*>& rods)
0028     : RodBarrelLayer(false), theRods(rods), theBinFinder(nullptr), isOverlapping(false) {
0029   // Sort rods in phi
0030   precomputed_value_sort(theRods.begin(), theRods.end(), geomsort::ExtractPhi<DetRod, float>());
0031
0032   theComponents.reserve(theRods.size());
0033   std::copy(theRods.begin(), theRods.end(), back_inserter(theComponents));
0034
0035   // Cache chamber pointers (the basic components_)
0036   for (const auto& it : rods) {
0037     vector<const GeomDet*> tmp2 = it->basicComponents();
0038     theBasicComps.insert(theBasicComps.end(), tmp2.begin(), tmp2.end());
0039   }
0040
0041   // Initialize the binfinder
0042   PhiBorderFinder bf(theRods);
0043   isOverlapping = bf.isPhiOverlapping();
0044
0045   if (bf.isPhiPeriodic()) {
0046     theBinFinder = new PeriodicBinFinderInPhi<double>(theRods.front()->position().phi(), theRods.size());
0047   } else {
0048     theBinFinder = new GeneralBinFinderInPhi<double>(bf);
0049   }
0050
0051   // Compute the layer's surface and bounds (from the components())
0052   BarrelDetLayer::initialize();
0053
0054   LogTrace("MTDDetLayers") << "Constructing MTDTrayBarrelLayer: " << basicComponents().size() << " Dets "
0055                            << theRods.size() << " Rods "
0056                            << " R: " << specificSurface().radius() << " Per.: " << bf.isPhiPeriodic()
0057                            << " Overl.: " << isOverlapping;
0058 }
0059
0060 MTDTrayBarrelLayer::~MTDTrayBarrelLayer() {
0061   delete theBinFinder;
0062   for (auto& i : theRods) {
0063     delete i;
0064   }
0065 }
0066
0067 vector<GeometricSearchDet::DetWithState> MTDTrayBarrelLayer::compatibleDets(
0068     const TrajectoryStateOnSurface& startingState, const Propagator& prop, const MeasurementEstimator& est) const {
0069   vector<DetWithState> result;
0070
0071   LogTrace("MTDDetLayers") << "MTDTrayBarrelLayer::compatibleDets, Cyl R: " << specificSurface().radius()
0072                            << " TSOS at R= " << startingState.globalPosition().perp()
0073                            << " phi= " << startingState.globalPosition().phi();
0074
0075   pair<bool, TrajectoryStateOnSurface> compat = compatible(startingState, prop, est);
0076   if (!compat.first) {
0077     LogTrace("MTDDetLayers") << "     MTDTrayBarrelLayer::compatibleDets: not compatible"
0078                              << " (should not have been selected!)";
0079     return vector<DetWithState>();
0080   }
0081
0082   TrajectoryStateOnSurface& tsos = compat.second;
0083
0084   LogTrace("MTDDetLayers") << "     MTDTrayBarrelLayer::compatibleDets, reached layer at: " << tsos.globalPosition()
0085                            << " R = " << tsos.globalPosition().perp() << " phi = " << tsos.globalPosition().phi();
0086
0087   int closest = theBinFinder->binIndex(tsos.globalPosition().phi());
0088   const DetRod* closestRod = theRods[closest];
0089
0090   // Check the closest rod
0091   LogTrace("MTDDetLayers") << "     MTDTrayBarrelLayer::compatibleDets, closestRod: " << closest
0092                            << " phi : " << closestRod->surface().position().phi()
0093                            << " FTS phi: " << tsos.globalPosition().phi();
0094
0095   result = closestRod->compatibleDets(tsos, prop, est);
0096
0097 #ifdef EDM_ML_DEBUG
0098   int nclosest = result.size();  // Debug counter
0099 #endif
0100
0101   bool checknext = false;
0102   double dist;
0103
0104   if (!result.empty()) {
0105     // Check if the track go outside closest rod, then look for closest.
0106     TrajectoryStateOnSurface& predictedState = result.front().second;
0107     float xErr = xError(predictedState, est);
0108     float halfWid = closestRod->surface().bounds().width() / 2.;
0109     dist = predictedState.localPosition().x();
0110
0111     // If the layer is overlapping, additionally reduce halfWid by 10%
0112     // to account for overlap.
0113     // FIXME: should we account for the real amount of overlap?
0114     if (isOverlapping)
0115       halfWid *= 0.9;
0116
0117     if (fabs(dist) + xErr > halfWid) {
0118       checknext = true;
0119     }
0120   } else {  // Rod is not compatible
0121     //FIXME: Usually next cannot be either. Implement proper logic.
0122     // (in general at least one rod should be when this method is called by
0123     // compatibleDets() which calls compatible())
0124     checknext = true;
0125
0126     // Look for the next-to closest in phi.
0127     // Note Geom::Phi, subtraction is pi-border-safe
0128     if (tsos.globalPosition().phi() - closestRod->surface().position().phi() > 0.) {
0129       dist = -1.;
0130     } else {
0131       dist = +1.;
0132     }
0133
0134     LogTrace("MTDDetLayers") << "     MTDTrayBarrelLayer::fastCompatibleDets, none on closest rod!";
0135   }
0136
0137   if (checknext) {
0138     int next;
0139     if (dist < 0.)
0140       next = closest + 1;
0141     else
0142       next = closest - 1;
0143
0144     next = theBinFinder->binIndex(next);  // Bin Periodicity
0145     const DetRod* nextRod = theRods[next];
0146
0147     LogTrace("MTDDetLayers") << "     MTDTrayBarrelLayer::fastCompatibleDets, next-to closest"
0148                              << " rod: " << next << " dist " << dist << " phi : " << nextRod->surface().position().phi()
0149                              << " FTS phi: " << tsos.globalPosition().phi();
0150
0151     vector<DetWithState> nextRodDets = nextRod->compatibleDets(tsos, prop, est);
0152     result.insert(result.end(), nextRodDets.begin(), nextRodDets.end());
0153   }
0154
0155 #ifdef EDM_ML_DEBUG
0156   LogTrace("MTDDetLayers") << "     MTDTrayBarrelLayer::fastCompatibleDets: found: " << result.size()
0157                            << " on closest: " << nclosest << " # checked rods: " << 1 + int(checknext);
0158 #endif
0159
0160   return result;
0161 }
0162
0163 vector<DetGroup> MTDTrayBarrelLayer::groupedCompatibleDets(const TrajectoryStateOnSurface& startingState,
0164                                                            const Propagator& prop,
0165                                                            const MeasurementEstimator& est) const {
0166   vector<GeometricSearchDet::DetWithState> detWithStates;
0167
0168   detWithStates = compatibleDets(startingState, prop, est);
0169
0170   vector<DetGroup> result;
0171   if (!detWithStates.empty()) {
0172     result.push_back(DetGroup(detWithStates));
0173   }
0174   LogTrace("MTDDetLayers") << "MTDTrayBarrelLayer Compatible rods: " << result.size();
0175   return result;
0176 }
0177
0178 GeomDetEnumerators::SubDetector MTDTrayBarrelLayer::subDetector() const { return theBasicComps.front()->subDetector(); }
0179
0180 const vector<const GeometricSearchDet*>& MTDTrayBarrelLayer::components() const { return theComponents; }
0181
0182 float MTDTrayBarrelLayer::xError(const TrajectoryStateOnSurface& tsos, const MeasurementEstimator& est) const {
0183   const float nSigmas = 3.f;
0184   if (tsos.hasError()) {
0185     return nSigmas * sqrt(tsos.localError().positionError().xx());
0186   } else
0187     return nSigmas * 0.5;
0188 }