DetLayers/src/MuRodBarrelLayer.cc

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