TkDetLayers/src/PixelBlade.cc

0001 #include "PixelBlade.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "BladeShapeBuilderFromDet.h"
0006 #include "LayerCrossingSide.h"
0007 #include "DetGroupMerger.h"
0008 #include "CompatibleDetToGroupAdder.h"
0009
0010 #include "TrackingTools/DetLayers/interface/DetLayerException.h"
0011 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0012 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
0013
0014 using namespace std;
0015
0016 typedef GeometricSearchDet::DetWithState DetWithState;
0017
0018 PixelBlade::~PixelBlade() {}
0019
0020 PixelBlade::PixelBlade(vector<const GeomDet*>& frontDets, vector<const GeomDet*>& backDets)
0021     : GeometricSearchDet(true), theFrontDets(frontDets), theBackDets(backDets) {
0022   theDets.assign(theFrontDets.begin(), theFrontDets.end());
0023   theDets.insert(theDets.end(), theBackDets.begin(), theBackDets.end());
0024
0025   theDiskSector = BladeShapeBuilderFromDet::build(theDets);
0026   theFrontDiskSector = BladeShapeBuilderFromDet::build(theFrontDets);
0027   theBackDiskSector = BladeShapeBuilderFromDet::build(theBackDets);
0028
0029   //--------- DEBUG INFO --------------
0030   LogDebug("TkDetLayers") << "DEBUG INFO for PixelBlade";
0031   LogDebug("TkDetLayers") << "Blade z, perp, innerRadius, outerR: " << this->position().z() << " , "
0032                           << this->position().perp() << " , " << theDiskSector->innerRadius() << " , "
0033                           << theDiskSector->outerRadius();
0034
0035   for (vector<const GeomDet*>::const_iterator it = theFrontDets.begin(); it != theFrontDets.end(); it++) {
0036     LogDebug("TkDetLayers") << "frontDet phi,z,r: " << (*it)->position().phi() << " , " << (*it)->position().z()
0037                             << " , " << (*it)->position().perp();
0038     ;
0039   }
0040
0041   for (vector<const GeomDet*>::const_iterator it = theBackDets.begin(); it != theBackDets.end(); it++) {
0042     LogDebug("TkDetLayers") << "backDet phi,z,r: " << (*it)->position().phi() << " , " << (*it)->position().z() << " , "
0043                             << (*it)->position().perp();
0044   }
0045   //-----------------------------------
0046 }
0047
0048 const vector<const GeometricSearchDet*>& PixelBlade::components() const {
0049   throw DetLayerException("TOBRod doesn't have GeometricSearchDet components");
0050 }
0051
0052 pair<bool, TrajectoryStateOnSurface> PixelBlade::compatible(const TrajectoryStateOnSurface& ts,
0053                                                             const Propagator&,
0054                                                             const MeasurementEstimator&) const {
0055   edm::LogError("TkDetLayers") << "temporary dummy implementation of PixelBlade::compatible()!!";
0056   return pair<bool, TrajectoryStateOnSurface>();
0057 }
0058
0059 void PixelBlade::groupedCompatibleDetsV(const TrajectoryStateOnSurface& tsos,
0060                                         const Propagator& prop,
0061                                         const MeasurementEstimator& est,
0062                                         std::vector<DetGroup>& result) const {
0063   SubLayerCrossings crossings;
0064   crossings = computeCrossings(tsos, prop.propagationDirection());
0065   if (!crossings.isValid())
0066     return;
0067
0068   vector<DetGroup> closestResult;
0069   addClosest(tsos, prop, est, crossings.closest(), closestResult);
0070
0071   if (closestResult.empty()) {
0072     vector<DetGroup> nextResult;
0073     addClosest(tsos, prop, est, crossings.other(), nextResult);
0074     if (nextResult.empty())
0075       return;
0076
0077     DetGroupElement nextGel(nextResult.front().front());
0078     int crossingSide = LayerCrossingSide().endcapSide(nextGel.trajectoryState(), prop);
0079
0080     DetGroupMerger::orderAndMergeTwoLevels(
0081         std::move(closestResult), std::move(nextResult), result, crossings.closestIndex(), crossingSide);
0082   } else {
0083     DetGroupElement closestGel(closestResult.front().front());
0084     float window = computeWindowSize(closestGel.det(), closestGel.trajectoryState(), est);
0085
0086     searchNeighbors(tsos, prop, est, crossings.closest(), window, closestResult, false);
0087
0088     vector<DetGroup> nextResult;
0089     searchNeighbors(tsos, prop, est, crossings.other(), window, nextResult, true);
0090
0091     int crossingSide = LayerCrossingSide().endcapSide(closestGel.trajectoryState(), prop);
0092     DetGroupMerger::orderAndMergeTwoLevels(
0093         std::move(closestResult), std::move(nextResult), result, crossings.closestIndex(), crossingSide);
0094   }
0095 }
0096
0097 SubLayerCrossings PixelBlade::computeCrossings(const TrajectoryStateOnSurface& startingState,
0098                                                PropagationDirection propDir) const {
0099   HelixPlaneCrossing::PositionType startPos(startingState.globalPosition());
0100   HelixPlaneCrossing::DirectionType startDir(startingState.globalMomentum());
0101   double rho(startingState.transverseCurvature());
0102
0103   HelixArbitraryPlaneCrossing crossing(startPos, startDir, rho, propDir);
0104
0105   pair<bool, double> innerPath = crossing.pathLength(*theFrontDiskSector);
0106   if (!innerPath.first)
0107     return SubLayerCrossings();
0108
0109   GlobalPoint gInnerPoint(crossing.position(innerPath.second));
0110   //Code for use of binfinder
0111   //int innerIndex = theInnerBinFinder.binIndex(gInnerPoint.perp());
0112   //float innerDist = std::abs( theInnerBinFinder.binPosition(innerIndex) - gInnerPoint.z());
0113   int innerIndex = findBin(gInnerPoint.perp(), 0);
0114   float innerDist = std::abs(findPosition(innerIndex, 0).perp() - gInnerPoint.perp());
0115   SubLayerCrossing innerSLC(0, innerIndex, gInnerPoint);
0116
0117   pair<bool, double> outerPath = crossing.pathLength(*theBackDiskSector);
0118   if (!outerPath.first)
0119     return SubLayerCrossings();
0120
0121   GlobalPoint gOuterPoint(crossing.position(outerPath.second));
0122   //Code for use of binfinder
0123   //int outerIndex = theOuterBinFinder.binIndex(gOuterPoint.perp());
0124   //float outerDist = std::abs( theOuterBinFinder.binPosition(outerIndex) - gOuterPoint.perp());
0125   int outerIndex = findBin(gOuterPoint.perp(), 1);
0126   float outerDist = std::abs(findPosition(outerIndex, 1).perp() - gOuterPoint.perp());
0127   SubLayerCrossing outerSLC(1, outerIndex, gOuterPoint);
0128
0129   if (innerDist < outerDist) {
0130     return SubLayerCrossings(innerSLC, outerSLC, 0);
0131   } else {
0132     return SubLayerCrossings(outerSLC, innerSLC, 1);
0133   }
0134 }
0135
0136 bool PixelBlade::addClosest(const TrajectoryStateOnSurface& tsos,
0137                             const Propagator& prop,
0138                             const MeasurementEstimator& est,
0139                             const SubLayerCrossing& crossing,
0140                             vector<DetGroup>& result) const {
0141   const vector<const GeomDet*>& sBlade(subBlade(crossing.subLayerIndex()));
0142   return CompatibleDetToGroupAdder().add(*sBlade[crossing.closestDetIndex()], tsos, prop, est, result);
0143 }
0144
0145 float PixelBlade::computeWindowSize(const GeomDet* det,
0146                                     const TrajectoryStateOnSurface& tsos,
0147                                     const MeasurementEstimator& est) const {
0148   return est.maximalLocalDisplacement(tsos, det->surface()).x();
0149 }
0150
0151 void PixelBlade::searchNeighbors(const TrajectoryStateOnSurface& tsos,
0152                                  const Propagator& prop,
0153                                  const MeasurementEstimator& est,
0154                                  const SubLayerCrossing& crossing,
0155                                  float window,
0156                                  vector<DetGroup>& result,
0157                                  bool checkClosest) const {
0158   const GlobalPoint& gCrossingPos = crossing.position();
0159
0160   const vector<const GeomDet*>& sBlade(subBlade(crossing.subLayerIndex()));
0161
0162   int closestIndex = crossing.closestDetIndex();
0163   int negStartIndex = closestIndex - 1;
0164   int posStartIndex = closestIndex + 1;
0165
0166   if (checkClosest) {  // must decide if the closest is on the neg or pos side
0167     if (gCrossingPos.perp() < sBlade[closestIndex]->surface().position().perp()) {
0168       posStartIndex = closestIndex;
0169     } else {
0170       negStartIndex = closestIndex;
0171     }
0172   }
0173
0174   typedef CompatibleDetToGroupAdder Adder;
0175   for (int idet = negStartIndex; idet >= 0; idet--) {
0176     if (!overlap(gCrossingPos, *sBlade[idet], window))
0177       break;
0178     if (!Adder::add(*sBlade[idet], tsos, prop, est, result))
0179       break;
0180   }
0181   for (int idet = posStartIndex; idet < static_cast<int>(sBlade.size()); idet++) {
0182     if (!overlap(gCrossingPos, *sBlade[idet], window))
0183       break;
0184     if (!Adder::add(*sBlade[idet], tsos, prop, est, result))
0185       break;
0186   }
0187 }
0188
0189 bool PixelBlade::overlap(const GlobalPoint& crossPoint, const GeomDet& det, float window) const {
0190   // check if the z window around TSOS overlaps with the detector theDet (with a 1% margin added)
0191
0192   //   const float tolerance = 0.1;
0193   const float relativeMargin = 1.01;
0194
0195   LocalPoint localCrossPoint(det.surface().toLocal(crossPoint));
0196   //   if (std::abs(localCrossPoint.z()) > tolerance) {
0197   //     edm::LogInfo(TkDetLayers) << "PixelBlade::overlap calculation assumes point on surface, but it is off by "
0198   //     << localCrossPoint.z() ;
0199   //   }
0200
0201   float localX = localCrossPoint.x();
0202   float detHalfLength = det.surface().bounds().length() / 2.;
0203
0204   //   edm::LogInfo(TkDetLayers) << "PixelBlade::overlap: Det at " << det.position() << " hit at " << localY
0205   //        << " Window " << window << " halflength "  << detHalfLength ;
0206
0207   if ((std::abs(localX) - window) < relativeMargin * detHalfLength) {  // FIXME: margin hard-wired!
0208     return true;
0209   } else {
0210     return false;
0211   }
0212 }
0213
0214 int PixelBlade::findBin(float R, int diskSectorIndex) const {
0215   vector<const GeomDet*> localDets = diskSectorIndex == 0 ? theFrontDets : theBackDets;
0216
0217   int theBin = 0;
0218   float rDiff = std::abs(R - localDets.front()->surface().position().perp());
0219   ;
0220   for (vector<const GeomDet*>::const_iterator i = localDets.begin(); i != localDets.end(); i++) {
0221     float testDiff = std::abs(R - (**i).surface().position().perp());
0222     if (testDiff < rDiff) {
0223       rDiff = testDiff;
0224       theBin = i - localDets.begin();
0225     }
0226   }
0227   return theBin;
0228 }
0229
0230 GlobalPoint PixelBlade::findPosition(int index, int diskSectorType) const {
0231   vector<const GeomDet*> diskSector = diskSectorType == 0 ? theFrontDets : theBackDets;
0232   return (diskSector[index])->surface().position();
0233 }