DetLayers/src/MuDetRing.cc

0001 /** \file
0002  *
0003  *  \author N. Amapane - CERN
0004  */
0005
0006 #include "RecoMuon/DetLayers/interface/MuDetRing.h"
0007 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0008 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0009 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 #include <iostream>
0013 #include <vector>
0014
0015 using namespace std;
0016
0017 MuDetRing::MuDetRing(vector<const GeomDet*>::const_iterator first, vector<const GeomDet*>::const_iterator last)
0018     : ForwardDetRingOneZ(first, last) {
0019   init();
0020 }
0021
0022 MuDetRing::MuDetRing(const vector<const GeomDet*>& vdets) : ForwardDetRingOneZ(vdets) { init(); }
0023
0024 void MuDetRing::init() {
0025   theBinFinder = BinFinderType(basicComponents().front()->position().phi(), basicComponents().size());
0026 }
0027
0028 MuDetRing::~MuDetRing() {}
0029
0030 const vector<const GeometricSearchDet*>& MuDetRing::components() const {
0031   // FIXME dummy impl.
0032   cout << "temporary dummy implementation of MuDetRing::components()!!" << endl;
0033   static const vector<const GeometricSearchDet*> result;
0034   return result;
0035 }
0036
0037 pair<bool, TrajectoryStateOnSurface> MuDetRing::compatible(const TrajectoryStateOnSurface& ts,
0038                                                            const Propagator& prop,
0039                                                            const MeasurementEstimator& est) const {
0040   const std::string metname = "Muon|RecoMuon|RecoMuonDetLayers|MuDetRing";
0041   TrajectoryStateOnSurface ms = prop.propagate(ts, specificSurface());
0042
0043   LogTrace(metname) << "MuDetRing::compatible, Surface at Z: " << specificSurface().position().z()
0044                     << " R1: " << specificSurface().innerRadius() << " R2: " << specificSurface().outerRadius()
0045                     << " TS   at Z,R: " << ts.globalPosition().z() << "," << ts.globalPosition().perp();
0046   if (ms.isValid()) {
0047     LogTrace(metname) << " DEST at Z,R: " << ms.globalPosition().z() << "," << ms.globalPosition().perp()
0048                       << " local Z: " << ms.localPosition().z() << endl;
0049   } else
0050     LogTrace(metname) << " DEST: not valid" << endl;
0051
0052   if (ms.isValid())
0053     return make_pair(est.estimate(ms, specificSurface()) != 0, ms);
0054   else
0055     return make_pair(false, ms);
0056 }
0057
0058 vector<GeometricSearchDet::DetWithState> MuDetRing::compatibleDets(const TrajectoryStateOnSurface& startingState,
0059                                                                    const Propagator& prop,
0060                                                                    const MeasurementEstimator& est) const {
0061   const std::string metname = "Muon|RecoMuon|RecoMuonDetLayers|MuDetRing";
0062
0063   LogTrace(metname) << "MuDetRing::compatibleDets, Surface at Z: " << surface().position().z()
0064                     << " R1: " << specificSurface().innerRadius() << " R2: " << specificSurface().outerRadius()
0065                     << " TS at Z,R: " << startingState.globalPosition().z() << ","
0066                     << startingState.globalPosition().perp() << "     DetRing pos." << position();
0067
0068   vector<DetWithState> result;
0069
0070   // Propagate and check that the result is within bounds
0071   pair<bool, TrajectoryStateOnSurface> compat = compatible(startingState, prop, est);
0072   if (!compat.first) {
0073     LogTrace(metname) << "    MuDetRing::compatibleDets: not compatible"
0074                       << "    (should not have been selected!)";
0075     return result;
0076   }
0077
0078   // Find the most probable destination component
0079   TrajectoryStateOnSurface& tsos = compat.second;
0080   GlobalPoint startPos = tsos.globalPosition();
0081   int closest = theBinFinder.binIndex(startPos.phi());
0082   const vector<const GeomDet*> dets = basicComponents();
0083   LogTrace(metname) << "     MuDetRing::compatibleDets, closest det: " << closest
0084                     << " Phi: " << dets[closest]->surface().position().phi() << " impactPhi " << startPos.phi();
0085
0086   // Add this detector, if it is compatible
0087   // NOTE: add performs a null propagation
0088   add(closest, result, tsos, prop, est);
0089
0090   int nclosest = result.size();
0091   int nnextdet = 0;  // MDEBUG counters
0092
0093   // Try the neighbors on each side until no more compatible.
0094   float dphi = 0;
0095   if (!result.empty()) {  // If closest is not compatible the next cannot be either
0096     float nSigmas = 3.;
0097     if (result.back().second.hasError()) {
0098       dphi = nSigmas * atan(sqrt(result.back().second.localError().positionError().xx()) /
0099                             result.back().second.globalPosition().perp());
0100     }
0101   } else {
0102     LogTrace(metname) << "     MuDetRing::compatibleDets, closest not compatible!";
0103     //FIXME:  if closest is not compatible the next cannot be either
0104   }
0105
0106   for (int idet = closest + 1; idet < closest + int(dets.size()) / 4 + 1; idet++) {
0107     // FIXME: should use dphi to decide if det must be queried.
0108     // Right now query until not compatible.
0109     int idetp = theBinFinder.binIndex(idet);
0110     {
0111       LogTrace(metname) << "     next det:" << idetp << " at Z: " << dets[idetp]->position().z()
0112                         << " phi: " << dets[idetp]->position().phi() << " FTS phi " << startPos.phi() << " max dphi "
0113                         << dphi;
0114       nnextdet++;
0115       if (!add(idetp, result, tsos, prop, est))
0116         break;
0117     }
0118   }
0119
0120   for (int idet = closest - 1; idet > closest - int(dets.size()) / 4 - 1; idet--) {
0121     // FIXME: should use dphi to decide if det must be queried.
0122     // Right now query until not compatible.
0123     int idetp = theBinFinder.binIndex(idet);
0124     {
0125       LogTrace(metname) << "     previous det:" << idetp << " " << idet << " " << closest - dets.size() / 4 - 1
0126                         << " at Z: " << dets[idetp]->position().z() << " phi: " << dets[idetp]->position().phi()
0127                         << " FTS phi " << startPos.phi() << " max dphi" << dphi;
0128       nnextdet++;
0129       if (!add(idetp, result, tsos, prop, est))
0130         break;
0131     }
0132   }
0133
0134   LogTrace(metname) << "     MuDetRing::compatibleDets, size: " << result.size() << " on closest: " << nclosest
0135                     << " # checked dets: " << nnextdet + 1;
0136
0137   if (result.empty()) {
0138     LogTrace(metname) << "   ***Ring not compatible,should have been discarded before!!!";
0139   }
0140
0141   return result;
0142 }
0143
0144 vector<DetGroup> MuDetRing::groupedCompatibleDets(const TrajectoryStateOnSurface& startingState,
0145                                                   const Propagator& prop,
0146                                                   const MeasurementEstimator& est) const {
0147   // FIXME should be implemented to allow returning  overlapping chambers
0148   // as separate groups!
0149   cout << "dummy implementation of MuDetRod::groupedCompatibleDets()" << endl;
0150   vector<DetGroup> result;
0151   return result;
0152 }