DetLayers/src/MuRingForwardDoubleLayer.cc

0001 /** \file
0002  *
0003  *  \author Rick Wilkinson
0004  */
0005
0006 #include <RecoMuon/DetLayers/interface/MuRingForwardDoubleLayer.h>
0007 #include <RecoMuon/DetLayers/interface/MuDetRing.h>
0008 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0009 #include <DataFormats/GeometrySurface/interface/SimpleDiskBounds.h>
0010 #include <TrackingTools/GeomPropagators/interface/Propagator.h>
0011 #include <TrackingTools/DetLayers/interface/MeasurementEstimator.h>
0012
0013 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0014
0015 #include <algorithm>
0016 #include <iostream>
0017 #include <vector>
0018
0019 using namespace std;
0020
0021 MuRingForwardDoubleLayer::MuRingForwardDoubleLayer(const vector<const ForwardDetRing*>& frontRings,
0022                                                    const vector<const ForwardDetRing*>& backRings)
0023     : RingedForwardLayer(true),
0024       theFrontLayer(frontRings),
0025       theBackLayer(backRings),
0026       theRings(frontRings),  // add back later
0027       theComponents(),
0028       theBasicComponents() {
0029   const std::string metname = "Muon|RecoMuon|RecoMuonDetLayers|MuRingForwardDoubleLayer";
0030
0031   theRings.insert(theRings.end(), backRings.begin(), backRings.end());
0032   theComponents = std::vector<const GeometricSearchDet*>(theRings.begin(), theRings.end());
0033
0034   // Cache chamber pointers (the basic components_)
0035   // and find extension in R and Z
0036   for (vector<const ForwardDetRing*>::const_iterator it = theRings.begin(); it != theRings.end(); it++) {
0037     vector<const GeomDet*> tmp2 = (*it)->basicComponents();
0038     theBasicComponents.insert(theBasicComponents.end(), tmp2.begin(), tmp2.end());
0039   }
0040
0041   setSurface(computeSurface());
0042
0043   LogTrace(metname) << "Constructing MuRingForwardDoubleLayer: " << basicComponents().size() << " Dets "
0044                     << theRings.size() << " Rings "
0045                     << " Z: " << specificSurface().position().z() << " R1: " << specificSurface().innerRadius()
0046                     << " R2: " << specificSurface().outerRadius();
0047
0048   selfTest();
0049 }
0050
0051 BoundDisk* MuRingForwardDoubleLayer::computeSurface() {
0052   const BoundDisk& frontDisk = theFrontLayer.specificSurface();
0053   const BoundDisk& backDisk = theBackLayer.specificSurface();
0054
0055   float rmin = min(frontDisk.innerRadius(), backDisk.innerRadius());
0056   float rmax = max(frontDisk.outerRadius(), backDisk.outerRadius());
0057   float zmin = frontDisk.position().z();
0058   float halfThickness = frontDisk.bounds().thickness() / 2.;
0059   zmin = (zmin > 0) ? zmin - halfThickness : zmin + halfThickness;
0060   float zmax = backDisk.position().z();
0061   halfThickness = backDisk.bounds().thickness() / 2.;
0062   zmax = (zmax > 0) ? zmax + halfThickness : zmax - halfThickness;
0063   float zPos = (zmax + zmin) / 2.;
0064   PositionType pos(0., 0., zPos);
0065   RotationType rot;
0066
0067   return new BoundDisk(pos, rot, new SimpleDiskBounds(rmin, rmax, zmin - zPos, zmax - zPos));
0068 }
0069
0070 bool MuRingForwardDoubleLayer::isInsideOut(const TrajectoryStateOnSurface& tsos) const {
0071   return tsos.globalPosition().basicVector().dot(tsos.globalMomentum().basicVector()) > 0;
0072 }
0073
0074 std::pair<bool, TrajectoryStateOnSurface> MuRingForwardDoubleLayer::compatible(
0075     const TrajectoryStateOnSurface& startingState, const Propagator& prop, const MeasurementEstimator& est) const {
0076   // mostly copied from ForwardDetLayer, except propagates to closest surface,
0077   // not to center
0078   const std::string metname = "Muon|RecoMuon|RecoMuonDetLayers|MuRingForwardDoubleLayer";
0079
0080   bool insideOut = isInsideOut(startingState);
0081   const MuRingForwardLayer& closerLayer = (insideOut) ? theFrontLayer : theBackLayer;
0082   LogTrace("Muon|RecoMuon|RecoMuonDetLayers|MuRingForwardDoubleLayer")
0083       << "MuRingForwardDoubleLayer::compatible is assuming inside-out direction: " << insideOut;
0084
0085   //std::pair<bool, TrajectoryStateOnSurface> result
0086   //  = closerLayer.compatible(startingState, prop, est);
0087   //if(!result.first)
0088   // {
0089   //  result = furtherLayer.compatible(startingState, prop, est);
0090   //}
0091
0092   TrajectoryStateOnSurface myState = prop.propagate(startingState, closerLayer.specificSurface());
0093   if (!myState.isValid())
0094     return make_pair(false, myState);
0095
0096   // take into account the thickness of the layer
0097   float deltaR = surface().bounds().thickness() / 2. * fabs(tan(myState.localDirection().theta()));
0098
0099   // take into account the error on the predicted state
0100   const float nSigma = 3.;
0101   if (myState.hasError()) {
0102     LocalError err = myState.localError().positionError();
0103     // ignore correlation for the moment...
0104     deltaR += nSigma * sqrt(err.xx() + err.yy());
0105   }
0106
0107   float zPos = (zmax() + zmin()) / 2.;
0108   SimpleDiskBounds tmp(rmin() - deltaR, rmax() + deltaR, zmin() - zPos, zmax() - zPos);
0109
0110   return make_pair(tmp.inside(myState.localPosition()), myState);
0111 }
0112
0113 vector<GeometricSearchDet::DetWithState> MuRingForwardDoubleLayer::compatibleDets(
0114     const TrajectoryStateOnSurface& startingState, const Propagator& prop, const MeasurementEstimator& est) const {
0115   vector<DetWithState> result;
0116   const std::string metname = "Muon|RecoMuon|RecoMuonDetLayers|MuRingForwardDoubleLayer";
0117   pair<bool, TrajectoryStateOnSurface> compat = compatible(startingState, prop, est);
0118
0119   if (!compat.first) {
0120     LogTrace(metname) << "     MuRingForwardDoubleLayer::compatibleDets: not compatible"
0121                       << " (should not have been selected!)";
0122     return result;
0123   }
0124
0125   TrajectoryStateOnSurface& tsos = compat.second;
0126
0127   // standard implementation of compatibleDets() for class which have
0128   // groupedCompatibleDets implemented.
0129   // This code should be moved in a common place intead of being
0130   // copied many times.
0131   vector<DetGroup> vectorGroups = groupedCompatibleDets(tsos, prop, est);
0132   for (vector<DetGroup>::const_iterator itDG = vectorGroups.begin(); itDG != vectorGroups.end(); itDG++) {
0133     for (vector<DetGroupElement>::const_iterator itDGE = itDG->begin(); itDGE != itDG->end(); itDGE++) {
0134       result.push_back(DetWithState(itDGE->det(), itDGE->trajectoryState()));
0135     }
0136   }
0137   return result;
0138 }
0139
0140 vector<DetGroup> MuRingForwardDoubleLayer::groupedCompatibleDets(const TrajectoryStateOnSurface& startingState,
0141                                                                  const Propagator& prop,
0142                                                                  const MeasurementEstimator& est) const {
0143   const std::string metname = "Muon|RecoMuon|RecoMuonDetLayers|MuRingForwardDoubleLayer";
0144   vector<GeometricSearchDet::DetWithState> detWithStates1, detWithStates2;
0145
0146   LogTrace(metname) << "groupedCompatibleDets are currently given always in inside-out order";
0147   // this should be fixed either in RecoMuon/MeasurementDet/MuonDetLayerMeasurements or
0148   // RecoMuon/DetLayers/MuRingForwardDoubleLayer
0149   // and removed the reverse operation in StandAloneMuonFilter::findBestMeasurements
0150
0151   detWithStates1 = theFrontLayer.compatibleDets(startingState, prop, est);
0152   detWithStates2 = theBackLayer.compatibleDets(startingState, prop, est);
0153
0154   vector<DetGroup> result;
0155   if (!detWithStates1.empty())
0156     result.push_back(DetGroup(detWithStates1));
0157   if (!detWithStates2.empty())
0158     result.push_back(DetGroup(detWithStates2));
0159   LogTrace(metname) << "DoubleLayer Compatible dets: " << result.size();
0160   return result;
0161 }
0162
0163 bool MuRingForwardDoubleLayer::isCrack(const GlobalPoint& gp) const {
0164   const std::string metname = "Muon|RecoMuon|RecoMuonDetLayers|MuRingForwardDoubleLayer";
0165   // approximate
0166   bool result = false;
0167   double r = gp.perp();
0168   const std::vector<const ForwardDetRing*>& backRings = theBackLayer.rings();
0169   if (backRings.size() > 1) {
0170     const MuDetRing* innerRing = dynamic_cast<const MuDetRing*>(backRings[0]);
0171     const MuDetRing* outerRing = dynamic_cast<const MuDetRing*>(backRings[1]);
0172     assert(innerRing && outerRing);
0173     float crackInner = innerRing->specificSurface().outerRadius();
0174     float crackOuter = outerRing->specificSurface().innerRadius();
0175     LogTrace(metname) << "In a crack:" << crackInner << " " << r << " " << crackOuter;
0176     if (r > crackInner && r < crackOuter)
0177       return true;
0178   }
0179   // non-overlapping rings
0180   //double phi = gp.phi().degrees();
0181   return result;
0182 }
0183
0184 void MuRingForwardDoubleLayer::selfTest() const {
0185   const std::vector<const GeomDet*>& frontDets = theFrontLayer.basicComponents();
0186   const std::vector<const GeomDet*>& backDets = theBackLayer.basicComponents();
0187
0188   std::vector<const GeomDet*>::const_iterator frontItr = frontDets.begin(), lastFront = frontDets.end(),
0189                                               backItr = backDets.begin(), lastBack = backDets.end();
0190
0191   // test that each front z is less than each back z
0192   for (; frontItr != lastFront; ++frontItr) {
0193     float frontz = fabs((**frontItr).surface().position().z());
0194     for (; backItr != lastBack; ++backItr) {
0195       float backz = fabs((**backItr).surface().position().z());
0196       assert(frontz < backz);
0197     }
0198   }
0199 }