TkDetLayers/src/TIBLayer.cc

0001 #include "TIBLayer.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "DataFormats/GeometrySurface/interface/SimpleCylinderBounds.h"
0006
0007 #include "LayerCrossingSide.h"
0008 #include "DetGroupMerger.h"
0009 #include "CompatibleDetToGroupAdder.h"
0010
0011 #include "TrackingTools/DetLayers/interface/DetLayerException.h"
0012 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0013 #include "TrackingTools/GeomPropagators/interface/HelixBarrelCylinderCrossing.h"
0014 #include "TrackingTools/DetLayers/interface/DetLessZ.h"
0015
0016 typedef GeometricSearchDet::DetWithState DetWithState;
0017
0018 TIBLayer::TIBLayer(std::vector<const TIBRing*>& innerRings, std::vector<const TIBRing*>& outerRings)
0019     : TBLayer(innerRings, outerRings, GeomDetEnumerators::TIB) {
0020   theComps.assign(theInnerComps.begin(), theInnerComps.end());
0021   theComps.insert(theComps.end(), theOuterComps.begin(), theOuterComps.end());
0022
0023   std::sort(theComps.begin(), theComps.end(), isDetLessZ);
0024   std::sort(theInnerComps.begin(), theInnerComps.end(), isDetLessZ);
0025   std::sort(theOuterComps.begin(), theOuterComps.end(), isDetLessZ);
0026
0027   for (auto& it : theComps) {
0028     theBasicComps.insert(theBasicComps.end(), it->basicComponents().begin(), it->basicComponents().end());
0029   }
0030
0031   // initialize the surface
0032   theInnerCylinder = cylinder(theInnerComps);
0033   theOuterCylinder = cylinder(theOuterComps);
0034   initialize();
0035
0036   LogDebug("TkDetLayers") << "==== DEBUG TIBLayer =====";
0037   LogDebug("TkDetLayers") << "innerCyl radius, thickness, lenght: " << theInnerCylinder->radius() << " , "
0038                           << theInnerCylinder->bounds().thickness() << " , " << theInnerCylinder->bounds().length();
0039
0040   LogDebug("TkDetLayers") << "outerCyl radius, thickness, lenght: " << theOuterCylinder->radius() << " , "
0041                           << theOuterCylinder->bounds().thickness() << " , " << theOuterCylinder->bounds().length();
0042
0043   LogDebug("TkDetLayers") << "Cyl radius, thickness, lenght: " << specificSurface().radius() << " , "
0044                           << specificSurface().bounds().thickness() << " , " << specificSurface().bounds().length();
0045
0046   for (auto& i : theInnerComps) {
0047     LogDebug("TkDetLayers") << "inner TIBRing pos z,radius,eta,phi: " << i->position().z() << " , "
0048                             << i->position().perp() << " , " << i->position().eta() << " , " << i->position().phi();
0049   }
0050
0051   for (auto& i : theOuterComps) {
0052     LogDebug("TkDetLayers") << "outer TIBRing pos z,radius,eta,phi: " << i->position().z() << " , "
0053                             << i->position().perp() << " , " << i->position().eta() << " , " << i->position().phi();
0054   }
0055
0056   // initialise the bin finders
0057   //  vector<const GeometricSearchDet*> tmpIn;
0058   //for (vector<const TIBRing*>::const_iterator i=theInnerRings.begin();
0059   //     i != theInnerRings.end(); i++) tmpIn.push_back(*i);
0060   theInnerBinFinder = GeneralBinFinderInZforGeometricSearchDet<float>(theInnerComps.begin(), theInnerComps.end());
0061
0062   theOuterBinFinder = GeneralBinFinderInZforGeometricSearchDet<float>(theOuterComps.begin(), theOuterComps.end());
0063 }
0064
0065 TIBLayer::~TIBLayer() {}
0066
0067 BoundCylinder* TIBLayer::cylinder(const std::vector<const GeometricSearchDet*>& rings) {
0068   float leftPos = rings.front()->surface().position().z();
0069   float rightPos = rings.back()->surface().position().z();
0070
0071   const BoundCylinder& frontRing = static_cast<const BoundCylinder&>(rings.front()->surface());
0072   const BoundCylinder& backRing = static_cast<const BoundCylinder&>(rings.back()->surface());
0073   float r = frontRing.radius();
0074   const Bounds& leftBounds = frontRing.bounds();
0075   const Bounds& rightBounds = backRing.bounds();
0076
0077   //float r = rings.front()->specificSurface().radius();
0078   //const Bounds& leftBounds = rings.front()->specificSurface().bounds();
0079   //const Bounds& rightBounds = rings.back()->specificSurface().bounds();
0080
0081   float thick = leftBounds.thickness() / 2;
0082   float zmin = leftPos - leftBounds.length() / 2;
0083   float zmax = rightPos + rightBounds.length() / 2;
0084   float rmin = r - thick;
0085   float rmax = r + thick;
0086   float zpos = 0.5 * (leftPos + rightPos);
0087
0088   auto scp = new SimpleCylinderBounds(rmin, rmax, zmin - zpos, zmax - zpos);
0089   return new Cylinder(r, Surface::PositionType(0, 0, zpos), rings.front()->surface().rotation(), scp);
0090 }
0091
0092 std::tuple<bool, int, int> TIBLayer::computeIndexes(GlobalPoint gInnerPoint, GlobalPoint gOuterPoint) const {
0093   int innerIndex = theInnerBinFinder.binIndex(gInnerPoint.z());
0094   const GeometricSearchDet* innerRing = theInnerComps[innerIndex];
0095   float innerDist = std::abs(innerRing->surface().position().z() - gInnerPoint.z());
0096
0097   int outerIndex = theOuterBinFinder.binIndex(gOuterPoint.z());
0098   const GeometricSearchDet* outerRing = theOuterComps[outerIndex];
0099   float outerDist = std::abs(outerRing->surface().position().z() - gOuterPoint.z());
0100
0101   return std::make_tuple(innerDist < outerDist, innerIndex, outerIndex);
0102 }
0103
0104 void TIBLayer::searchNeighbors(const TrajectoryStateOnSurface& tsos,
0105                                const Propagator& prop,
0106                                const MeasurementEstimator& est,
0107                                const SubLayerCrossing& crossing,
0108                                float window,
0109                                std::vector<DetGroup>& result,
0110                                bool checkClosest) const {
0111   const GlobalPoint& gCrossingPos = crossing.position();
0112
0113   const std::vector<const GeometricSearchDet*>& sLayer(subLayer(crossing.subLayerIndex()));
0114
0115   int closestIndex = crossing.closestDetIndex();
0116   int negStartIndex = closestIndex - 1;
0117   int posStartIndex = closestIndex + 1;
0118
0119   if (checkClosest) {  // must decide if the closest is on the neg or pos side
0120     if (gCrossingPos.z() < sLayer[closestIndex]->surface().position().z()) {
0121       posStartIndex = closestIndex;
0122     } else {
0123       negStartIndex = closestIndex;
0124     }
0125   }
0126
0127   typedef CompatibleDetToGroupAdder Adder;
0128   for (int idet = negStartIndex; idet >= 0; idet--) {
0129     const GeometricSearchDet* neighborRing = sLayer[idet];
0130     if (!overlap(gCrossingPos, *neighborRing, window))
0131       break;
0132     if (!Adder::add(*neighborRing, tsos, prop, est, result))
0133       break;
0134   }
0135   for (int idet = posStartIndex; idet < static_cast<int>(sLayer.size()); idet++) {
0136     const GeometricSearchDet* neighborRing = sLayer[idet];
0137     if (!overlap(gCrossingPos, *neighborRing, window))
0138       break;
0139     if (!Adder::add(*neighborRing, tsos, prop, est, result))
0140       break;
0141   }
0142 }
0143
0144 bool TIBLayer::overlap(const GlobalPoint& crossPoint, const GeometricSearchDet& det, float window) {
0145   float halfLength = 0.5f * det.surface().bounds().length();
0146
0147   //   edm::LogInfo(TkDetLayers) << " TIBLayer: checking ring with z " << det.position().z();
0148
0149   return std::abs(crossPoint.z() - det.position().z()) < (halfLength + window);
0150 }
0151
0152 float TIBLayer::computeWindowSize(const GeomDet* det,
0153                                   const TrajectoryStateOnSurface& tsos,
0154                                   const MeasurementEstimator& est) const {
0155   // we assume the outer and inner rings have about same thickness...
0156
0157   //   edm::LogInfo(TkDetLayers) << "TIBLayer::computeWindowSize: Y axis of tangent plane is"
0158   //        << plane.toGlobal( LocalVector(0,1,0)) ;
0159
0160   MeasurementEstimator::Local2DVector localError(est.maximalLocalDisplacement(tsos, det->surface()));
0161   float yError = localError.y();
0162
0163   // float tanTheta = std::tan( tsos.globalMomentum().theta());
0164   auto gm = tsos.globalMomentum();
0165   auto cotanTheta = gm.z() / gm.perp();
0166   float thickCorrection = 0.5f * det->surface().bounds().thickness() * std::abs(cotanTheta);
0167
0168   // FIXME: correct this in case of wide phi window !
0169
0170   return yError + thickCorrection;
0171 }